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Abstract 

The  growing  concern  over  the  impacts  of  oil  spills  on  aquatic  environments  has 
led  to  the  development  of  many  computer  models  for  simulating  the  transport 
and  spreading  of  oil  slicks  in  surface  waters.  Almost  oil  of  these  models  were 
developed  for  coastal  environments .  A  few  river  models  exist.  These  models  only 
considered  the  movement  of  surface  oil  slicks.  In  this  study  a  two-layer  model 
ROSS2,  is  developed  for  simulating  oil  spills  in  rivers.  This  model  considers  the 
oil  in  the  river  to  consist  of  a  surface  slick  and  suspended  oil  droplets  entroined 
over  the  depth  of  the  flow.  The  oil  transformation  processes  considered  in  the 
model  include  advection,  mechanical  spreading,  turbulent  diffusion  and  mix¬ 
ing,  evaporation,  dissolution,  emulsification,  shoreline  deposition  and  sinking. 
The  model  can  be  used  for  simulating  instantaneous  or  continuous  spills  either 
on  or  under  the  water  surfoce  in  rivers  with  or  without  an  ice  cover.  The  model 
has  been  implemented  for  the  Ohio-Monongahelo-Allegheny  river  system  and 
the  upper  St.  Lawrence  River.  This  report  describes  the  model  formulation  and 
implementation.  A  case  study  is  presented  along  with  detoiied  explanations  of 
the  program  structure  and  its  input  and  output.  Although  it  is  developed  for 
simulating  oil  spills,  the  model  can  be  applied  to  spills  of  other  hazardous 
materials. 
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A  Mathematical  Model  for 
Oil  Slick  Transport  and  Mixing  in  Rivers 

HUNG  TAG  SHEN,  POOJITHA  D.  YAPA,  DE  SHENG  WANG  AND  XIAO  QING  YANG 

INTRODUCTION 

In  recent  years  there  has  been  a  growing  concern  over  the  increasing  contamination  of  waterways 
and  shoreline  areas  caused  by  oil  spills.  CHI  spills  in  inland  waterways  can  have  enormous  environmen¬ 
tal  and  economic  impacts.  Oil  pollution  in  rivers  can  not  only  cause  long-tenn  damage  to  the  aquahc 
environment  for  fish  and  wildlife  but  also  pose  threats  to  water  supplies  for  areas  along  the  river.  With 
heavy  industrial  development  in  areas  along  rivers  and  inland  navigation  activities,  major  rivers  in  the 
United  States  are  vulnerable  to  oil  spills.  For  example,  in  January  1988,  an  oil  storage  tank  at  West 
Elizabeth,  Pennsylvania,  collapsed,  spilling  over  17,000  barrels  of  diesel  oil  into  the  Monongahela  Riv¬ 
er.  The  oil  sbck  reached  the  Ohio  River  near  Pittsburgh,  which  is  25  miles  downstream  of  the  spill  site, 
within  a  day  and  drifted  farther  downstream,  affecting  water  supply  intakes  along  the  river  (Miklaucic 
and  Saseen  1989). 

C(»icem  over  the  adverse  impacts  of  oil  spills  has  led  government  agencies  and  private  industries  to 
develop  oil  spill  emergency  response  plans  (e.g.  Hung  1991).  An  important  element  in  diese  programs 
is  the  use  of  computer  models  to  predict  the  movement  and  possible  impact  of  an  oil  spill.  In  the  event 
of  a  real  oil  ^ill,  a  model  can  be  used  on  a  real-time  basis  to  assist  dre  containment  and  recovery  of  the 
oil,  and  it  can  also  help  to  guide  the  field  data  collection  for  detailed  environmental  impact  analysis. 
These  objectives  can  be  achieved  by  using  die  model  to  forecast  the  location  and  distribution  of  the  oil. 
A  computer  model  can  also  be  used  to  study  scenarios  of  possible  spills  to  assist  in  developing  contin¬ 
gency  plans  and  assessing  likely  environmental  impacts. 

Numerous  oil  spill  models  have  been  developed  during  the  last  two  decades,  as  reviewed  by  Stolzen- 
badi  et  al.  (1977),  Huang  and  Monastero  (1982)  and  Spaulding  (1988).  Almost  all  of  these  models  were 
developed  for  coastal  marine  environments.  In  a  recent  study  Shen  and  Yapa  (1988)  developed  a  mod¬ 
el,  ROSS,  for  surfiice  oil  slkk  transport  in  rivers.  Oil  slick  transformation  processes  considered  in  the 
model  include  advection,  mechanical  spreading,  horizontal  turbulent  diffusion,  evaporation,  dissolu¬ 
tion  and  shoreline  deposition.  In  this  study,  a  two-layer  model,  ROSS2,  is  developed  to  take  into  con- 
irideration  tibe  di^>errion  of  oQ  droplets  into  the  water  column,  as  well  as  related  processes  including 
subsurface  tran^xirt,  mass  exchanges  between  the  surface  slick  and  the  suspended  oil  droplets  in  the 
water  column,  and  deposition  of  oil  to  he  channel  bottom  due  to  sedimentation. 

Tnoispoit  ai^  fate  of  spilled  oil 

In  addition  to  the  location,  size  and  physical-chemical  properties  of  he  spill,  he  fate  of  an  oil  spill  is 
infected  by  complac  intraielated  transport  and  weatitering  processes  (Mackay  and  McAuliffe  1988). 


Sinking 


^degrad^n 


Bottom  Sediment 

1  i 

1 

1 

> 

Hours  Days 

Weeks 

Months 

Years'" 

Time  From  Onset  of  Spill 

Figure  1.  Physical,  chemical  and  biological  processes  affecting  the  oil  slick  transformation. 


Figure  1  shows  a  schematic  representation  of  the  fate  of  oil  spilled  in  a  surface  water  body.  Oil  may 
enter  a  water  body  either  through  a  surface  spill  or  an  underwater  discharge.  A  part  of  the  spilled  oil 
will  form  a  surface  slick  on  the  water  surface  and  can  be  moved  about  by  the  action  of  winds,  currents 
or  waves.  Some  low-molecular-weight  hydrocarbons  can  dissolve  into  the  underlying  water  column, 
but  most  of  these  are  lost  to  the  atmosphere  through  evaporation.  Evaporation  of  volatile  components 
will  reduce  the  volume  of  spilled  oil  by  as  much  as  50%  during  the  first  few  days  after  the  spill  (Butler 
et  al.  1976).  In  turbulent  waters  some  of  the  oil  is  emulsified  and  dispersed  into  the  water  column  as 
suspended  droplets.  Some  of  the  oil  droplets  in  the  suspension  may  become  attached  to  suspended 
particulate  matter  and  slowly  settle  to  the  bottom;  some  will  rise  to  the  water  surface  due  to  buoyancy 
to  form  a  water-in-oU  emulsion.  Current  and  waves  may  also  drive  oil  onto  beaches  or  riverbanks.  For 
spills  in  rivers  the  transport  of  oil  is  also  affected  by  hydraulic  structures,  such  as  locks  and  dams,  as 
shown  in  Figure  2.  In  addition,  photochemical  reactions  and  microbial  biodegradation  can  change  the 
character  of  the  oil  and  reduce  the  amount  of  oil  present. 


Figure  2.  OU  slick  tranfforrmtion  in  rivers. 
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Advection 

Advection  is  a  physical  process  that  involves  drifting  of  the  surface  oil  slick  and  the  subsurface  oil.  It 
is  the  main  mechanism  that  governs  the  locarion  of  the  oil  following  its  discharge  into  the  river.  Once 
spilled  into  the  river,  the  oil  will  be  transported  downstream  by  the  flow,  with  the  influence  of  wind. 
The  advection  of  surface  oil  is  caused  by  the  combined  effects  of  surface  current  and  wind  drag.  The 
advection  of  subsurface  oil  is  the  movement  of  suspended  oil  droplets  entrained  in  the  flow  due  to  the 
subsurface  current. 

Spreading 

Spreading  is  the  physical  process  that  involves  drifting  of  the  surface  oil  slick  and  the  subsurface  oil 
as  it  is  being  transported  in  the  river.  The  spreading  phenomena  include  mechanical  spreading  and 
horizontal  diffusion.  Mechanical  spreading  is  the  horizontal  spreading  of  the  surface  oil  shck  due  to  the 
balancing  forces  of  inertia,  gravity,  viscosity  and  surface  tension.  This  mechanism  terminates  when  the 
thickness  of  the  slick  is  reduced  and  the  thin  oil  slick  ruptures  into  small  patches.  Horizontal  diffusion 
is  the  spreading  of  oil  due  to  the  turbulent  fluctuations  of  wind  and  current  velocities.  This  is  the  main 
mechanism  that  causes  the  spreading  of  both  the  surface  and  subsurface  oils.  Since  the  spreading  of  oil 
enhances  other  weather  processes,  such  as  evaporation,  dissolution  and  emulsification,  it  is  one  of  the 
most  important  processes  affecting  the  fate  of  the  spilled  oil. 

Evaporation 

Evaporation  occurs  immediately  after  the  spill.  As  the  surface  slick  spreads,  more  of  the  hydrocar¬ 
bons  are  exposed  to  the  atmosphere,  causing  the  evaporation  rate  to  increase.  The  amount  and  rate  of 
evaporation  depend  on  the  percentage  of  light,  or  volatile,  components  in  oil.  Evaporation  is  the  most 
significant  physical-chemical  process  causing  the  reduction  in  oil  volume.  Highly  refined  oil  can  lose 
75%  or  more  of  its  volume  through  evaporation  within  a  matter  of  days. 

Dissolution 

Some  of  the  same  oil  components  that  are  subjected  to  evaporation  can  also  dissolve  into  the  water 
column  from  a  surface  slick.  Only  the  low-molecular-weight  hydrocarbons  have  an  appreciable  solu¬ 
bility  in  water.  The  fraction  of  oil  dissolved  is  very  small  compared  to  the  oil  that  is  evaporated.  How¬ 
ever,  this  extraction  process  can  be  important  because  of  the  toxidty  of  the  dissolved  fraction. 

Emulsification 

Dispersion,  or  oil-in-water  emulsification,  is  the  mixing  of  surface  oil  into  the  water  column.  Break¬ 
ing  waves  and  other  surface  turbulence  can  mix  the  surface  oil  into  flie  water  column  by  forming  many 
anall  globules  of  oil  that  can  be  rapidly  dispersed  over  the  depth  of  the  river  flow  and  transported  by 
the  current.  Some  of  tite  fine  oil  droplets  in  the  water  column  may  dissolve  in  the  water  column,  attach 
to  solid  particles  and  biodegrade.  A  significant  portion  of  suspended  oil  droplets  will  rise  to  flte  water 
surface  when  the  buoyancy  force  is  large  enough  to  overcome  the  vertical  mixing.  This  resurfacing 
mechanism  can  increase  the  area  of  oil  on  the  water  surface  and  fire  duration  of  oil  passage  at  any 
particular  site  along  the  river. 

Water-in-oii  emulsions  can  also  be  formed,  particularly  with  heavy  crudes  and  residual  oil.  The  re¬ 
sulting  emulsicms  contain  a  large  percentage  of  water  but  have  a  semisolid  texture,  often  referred  to  as 
"mousse"  because  of  their  appearance.  The  formation  of  mousse  is  a  process  of  dispersion  of  emulsi¬ 
fied  water  into  the  oil.  This  process  will  increase  the  viscosity  and  volume  of  the  slick.  Stable  emulsions 
am  have  a  water  ccmtent  of  up  to  80%.  The  mechanism  of  water-in-oil  emulsification  is  not  clear,  al¬ 
though  it  is  believed  that  this  is  due  to  the  presence  of  asphaltertes,  waxes  and  surfactants  and  the 
turtodenoe  associated  with  dispersion  and  reWfadng  processes. 
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Photo-oxidation 

In  the  presence  of  atmospheric  oxygen,  natural  sunlight  has  sufficient  energy  to  change  the  compo¬ 
sition  of  the  oil.  This  photo-oxidation  process  is  very  slow  and  usually  becomes  important  a  few  days 
after  the  spill.  The  extent  and  rate  of  photo-oxidation  depends  primarily  on  the  chemical  composition 
and  optical  density  of  the  oil.  The  photochemical  reactions  can  change  the  interfacial  properties  of  the 
oil  and  may  dissolve  toxic  organic  species  into  the  water  column.  However,  its  effect  on  the  oil  slick 
transport  process  is  not  very  important. 

Shoreline  deposition  and  sedimentation 

Oil  may  be  brought  to  the  riverfaanks  and  deposited  along  the  shoreline,  to  be  later  re-entrained  into 
the  river  current.  This  process  can  significantly  affect  the  distribution  of  oil  and  should  be  modeled 
Some  of  the  suspended  oil  droplets  may  also  sink  to  the  riverbed.  This  sinking  or  sedimentation  pro¬ 
cess  occurs  due  to  an  increase  in  density  of  the  oil,  resulting  from  either  the  evaporation  and  dissolu¬ 
tion  of  lighter  fractiorrs  of  the  oil,  or  adherence  of  the  oil  droplets  onto  suspended  sediment.  The  oil 
deposited  on  the  channel  bottom  may  be  moved  laterally  or  resuspended,  or  it  may  undergo  further 
biological  or  physical-chemical  reaction.  Little  is  known  about  the  ultimate  fate  of  the  sedimented  oil. 

Biodegradation 

All  of  the  processes  just  described,  except  possibly  photo-oxidation,  can  oidy  redistribute  the  oil. 
They  cannot  remove  the  hydrocarbon  from  the  environment.  Real  degradation  takes  place  only  through 
biochemical  oxidation.  This  biodegradation  process,  which  may  continue  for  years  after  a  spill,  is  the 
principal  long-term  means  of  removing  the  spilled  oil  from  the  environment. 

Oil  spill  models 

Most  of  the  existing  oil  spill  models  simulate  only  the  advection  and  spreading  of  a  surface  slick. 
Other  models  deal  extensively  with  physical-chemical  processes  but  lack  the  component  for  simulat¬ 
ing  the  movement  of  the  slick.  Only  in  recent  models  have  the  incorporation  of  both  transport  and 
weathering  processes  been  attempted  (Huang  and  Monastero  1982,  Spaulding  1988).  Since  there  is  a 
significant  lack  of  data  for  a  reliable  anal)dical  formulation  to  be  established  for  many  of  the  weathering 
processes,  it  is  impractical  to  include  all  of  them  into  an  oil  spill  simulation  model.  It  would  be  more 
useful  to  include  the  most  significant  processes,  i.e.  those  accounting  for  the  bulk  of  the  oil,  while  omit¬ 
ting  others  so  that  uncertainty  in  the  outcome  can  be  reduced.  Almost  all  of  the  existing  oil  slick  models 
were  developed  for  coastal  marine  environments. 

Only  a  few  models  were  developed  for  inland  waters  (Huang  and  Monastero  1982).  Tsahalis  (1979b) 
developed  a  simulation  model  for  predicting  the  transport,  spreading  and  associated  shoreline  con¬ 
tamination  of  a  surface  oil  sUck  in  a  river.  In  his  model  the  current  velocity  distribution  in  the  river  is 
calculated  by  empirical  relationships  determined  from  field  data  with  some  modifications  for  the  sec¬ 
ondary  current  in  river  bends.  The  oil  slick  is  assumed  to  remam  circular,  wifii  its  radius  calculated 
according  to  Fa)r's  spreading  laws  (Fay  1969, 1971).  The  drift  velocity  of  the  slick  center  is  determined 
by  formulas  derived  by  Tsahalis  (1979a)  from  laboratory  experiments.  Fingas  and  Sydor  (1980)  devel¬ 
oped  a  two-dimensional  model  for  surface  oil  slicks  in  a  short  river  reach.  The  current  distribution  was 
determined  using  the  fiiute-difference  sdwme  of  Leendertse  (1970).  The  entire  oil  slick  voliune  is  repre¬ 
sented  by  a  large  number  of  individual  parcels.  The  drift  velocities  of  these  parcels  are  determined  by 
tire  wind  factca  approach.  A  random  fluctuation  component  is  included  to  represent  horizontal  diffu¬ 
sion.  The  area  of  each  oil  parcel  is  calculated  by  Fa/s  spreading  laws  for  circular  slicks. 

Recontly  9ien  and  Yapa  (1988)  developed  a  computer  model,  ROSS,  for  svirface  slicks  in  rivers.  The 
model  considered  flw  ef^  of  ice  covers.  In  this  m^el  the  oil  slick  is  considered  to  be  composed  of  a 
large  number  of  discrete  parcels  with  eciual  volume.  These  parcels  are  tracked  for  their  positions  and 
subjected  to  esn  e<pial  mass  loss  rate  during  the  simulation.  The  twodimensional  velocity  distribution 
of  flie  river  water  is  ccnnputed  using  a  separate  hydrodynamic  model.  The  advection  of  each  parcel  in 
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INPUT 


•  Model  parameters;  Evaporation,  dissolution,  dispersion,  resurfacing,  sedime.itation 

•  Oil  spill  information:  Spill  Ic.  uon  and  duration,  volume  and  piuperties  of  oil 

•  River  and  weather  datn  r^.^ometry.  wind  velocity,  discharge  and  water  levels,  air 
temperature,  ice  con'  ,,on,  river  coordinate  system  transformation 

I 

DETERMINE  DEPTH  AND  VELOCITY  DISTRIBUTIONS 

*  Two-dimensional  depth  distribution 

•  Two-dimensional  depth-averaged  velocity  distribution 

ADVECTION  AND  HORIZONTAL  DIFFUSION 

•  Surface  layer 

•  Suspended  layer 

i 

EMULSIFICATION.  RESURFACING  AND  SEDIMENTATION 

•  Some  surface  oil  is  emulsified  into  water 

•  Some  suspended  oil  rises  to  surface  and  some  sinks  to  the  bottom 

i 

EVAPORATION  AND  DISSOLUTION 
{only  for  surface  layer) 

SHORELINE  DEPOSITION 
(only  for  surface  layer) 

i 

MECHANICAL  SPREADING 
(only  for  surface  layer) 

I 

OUTPUT 

•  Volume  ^  oil  parcels  and  distributions  in  surface  and  suspended  layers 

•  Amount  and  location  of  deposition  on  shorelines 

•  Amounts  of  evaporation,  dissohition,  emulsification,  resurfacing  and  sedimentation 


Figure  3.  Structure  of  the  simulation  model 


ttie  slide  is  detennined  from  the  wind  velodty  and  the  computed  current  velodty,  using  the  wind  factor 
approadu  The  spreading  of  the  slick  is  simulated  by  considering  both  mechanical  spreading  and  hori¬ 
zontal  turbulent  diffusiotr.  The  model  developed  by  Fay  (1969, 1971)  is  used  to  simulate  mechanical 
^reading,  while  a  random-walk  simulation  is  used  to  simulate  spreading  due  to  surface  turbulence. 
The  loss  of  oil  due  to  the  shorelirre  deposition  is  calculated  according  to  the  oil-retention  capability  of 
tihe  shorelines.  The  losses  of  oil  due  to  evaporation  and  dissolution  are  calculated  based  on  empirical 
&Hxnulatkn\s  that  consider  tiie  effects  of  slick  area,  wind  velocity,  temperature  and  oil  properties.  Ncaie 
oi  the  above  models  cemadered  the  motion  of  oil  in  suspension  and  the  exchange  between  the  siu^ce 
oil  and  tite  suspended  oil.  Sedimentation  of  oil  to  the  bottom  was  also  not  considered. 

&i  tire  pteaeat  study  fte  model  ROSS  is  generalized  by  considering  the  movement  of  oil  both  on  and 
tinder  tite  water  surf»3e.  The  present  model,  ROSS2,  is  primarily  designed  for  tite  transport  of  oil  spilled 
Bt  rivers;  however,  it  can  al^  be  used  for  spills  of  ofoer  hazardous  materiafe.  A  brief  outline  of  the 
structure  of  the  txuxlel  is  presented  in  Figure  3. 
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In  RC)SS2  the  spilled  oil  is  divided  into  a  surface  layer  and  a  suspended  layer,  with  continuous 
exchanges  of  oO  mass  between  the  two  layers.  In  the  simulation,  oil  in  both  layers  is  discretized  into  a 
large  number  of  parcels  that  are  tracked  for  their  positions  and  volumes  at  each  time  step  during  the 
simulation.  Besides  the  two-dimensional  velocity  distribution,  the  distribution  of  How  depth  is  also 
determined  for  simulating  transport  in  the  suspended  layer.  Because  of  the  exchange  between  the  sur¬ 
face  layer  and  the  suspended  layer,  the  number  of  oil  parcels,  as  well  as  the  volume  of  each  parcel,  will 
change  from  time  to  time.  The  advection  of  oil  parcels  in  the  suspended  layer  is  determined  by  the 
depth-averaged  current  velocity,  and  those  on  the  surface  are  determined  by  the  wind  factor  approach. 
The  model  developed  by  Fay  (1969,  1971)  is  used  to  simulate  mechanical  spreading,  ^ind  a  random- 
walk  simulation  is  used  to  account  for  horizontal  spreading  due  to  turbulent  diffusion.  The  loss  of  oil 
in  the  surface  layer  due  to  evaporation  and  dissolution  is  calculated  based  on  analytical  and  empirical 
formulations,  which  consider  the  effects  of  slick  area,  wind  velocity,  temperature  and  oil  properties. 
The  loss  of  oil  due  to  shoreline  deposition  is  calculated  according  to  the  oil-retention  capability  of  the 
shorelines.  The  amounts  of  vertical  dispersion  into  the  water  column,  res'jrfacing  and  deposition  to  the 
channel  bottom  are  calculated  based  on  the  depth-averaged  oil  concentration  in  the  suspended  layer. 
Weathering  processes  that  occur  long  cifter  the  onset  of  the  spill  are  not  well  understood  and  less  signif¬ 
icant,  so  they  are  not  considered  in  this  model.  This  is  also  justified  from  the  operational  point  of  view, 
since  the  main  objechve  of  the  model  is  for  short-term  tactical  forecasts  to  assist  in  designing  response 
measures  in  the  event  of  a  major  spill 

The  present  simulation  model  has  been  applied  to  the  upper  St.  Lawrence  River  and  the  Ohio- 
Monongahela-AUegheny  river  system.  This  report  will  use  the  Ohio-Monongahela-Allegheny  river 
system  as  an  example  to  illustrate  the  model  implementation  in  detail.  The  Ohio-Monongahela-Al¬ 
legheny  river  system  is  more  complex  and  difficult  to  implement,  so  it  seiTcs  as  a  better  example. 


MODEL  FORMULATION 
Governing  equations 

In  the  present  model,  oil  in  the  river  is  considered  to  consist  of  a  surface  slick  coating  on  the  v.  ater 
surface  and  a  mixed  layer  containing  suspended  oil  droplets  extending  over  the  depth  of  the  flow,  with 
continuous  exchange  between  the  two  layers.  The  thickness  of  the  surface  layer  is  assumed  to  be  negli¬ 
gible  in  comparison  to  the  suspended  layer.  The  equation  of  motion  for  oil  in  the  surface  layer  can  be 
written  as; 


+  ai  Vb  Cv|  z  =0  -  T  Cs  -  Ca  Se  +  Ms  [x,  y)  -  Dg  {x,  y)  (1) 

where  x,y,t  =  space  and  time  variables 

2  =  vertical  coordinate  measured  downward  from  the  water  surface 
Cg  =  local  volumetric  oil  concentration  in  the  surface  layer  per  unit  surface  area 
Cy  =  local  volumetric  concentration  of  oil  in  the  suspended  layer  per  unit  volume  of  water 
Ug,  Vg  =  components  of  surface  drift  velocity  in  x  and  y  directions,  respectively 
Dx,  Dy,  =  difrusion  coefficients  in  the  x  and  y  directions,  respectively 

CTj  =  coefficient  representing  the  probability  of  deposition  of  an  oil  droplet  reaching  the 
water  surface 

Vj,  =  buoyant  velocity  of  suspended  oil  globules 
7  =  coefficient  describing  the  rate  at  which  the  surface  oil  is  dispersed  into  the  water 
column 

Q  =  area  concentration  of  oil  in  the  surface  layer 
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Sg  =  rate  of  evaporation  and  dissolution  per  unit  area  of  the  surface  slick 
Mj  =  effect  on  the  distribution  of  surface  oil  by  mechanical  spreading 
D5  =  effect  on  the  distribution  of  surface  oil  by  shoreline  deposition. 


The  concentration  distribution  in  the  suspended  layer  can  be  described  by: 
3 (Cyh )  ^  I  nix 


+  7Cs-  Pi t'vl 


where  C^.  =  depth-averaged  volumetric  concentration  of  oil  in  the  suspended  layer 
h-  flow  depth 

u,v=  components  of  depth-averaged  river  current  in  .v  and  1/  directions,  respectively 
Pi  =  coefficient  defining  the  rate  of  net  oil  deposition  on  the  riverbed  per  unit  area. 


(2) 


Equation  2  can  be  simplified  by  introducing  the  mass  conservation  equation  of  water: 
dt  dx  dy 


In  Lagrangian  form  the  simplified  equation  for  the  suspended  layer  becomes; 

Tv  +  ICs-pCv 

h 


Dt  ~  h 


IhD^^ 

.-i- 

>1 

1 

Q 

-  tc 

dx 

\  3x1 

h 

where  ^  =  2^  +  u  —  +  v~; 
Dt  dt  dx  dy 


aCy  —  ajCvl^  ; 


(3) 


(4) 


The  movement  of  spilled  oil  in  the  river  as  described  by  eq  1  and  4  is  mainly  governed  by  advechon 
and  diffusion.  These  equations  can  be  solved  when  the  flow  condition  is  known.  Finite-difference  and 
finite-element  methods  for  solving  advection-diffusion  equations  exist.  These  methods,  however,  of¬ 
ten  suffer  from  numerical  diffusion  problems  and  aemand  e.vcessi  ve  computational  time  and  memory' 
when  applied  to  long  river  reaches.  In  this  study  a  modified  Lagrangian  discrete-parcel  method  is  used. 
The  method  is  a  generalization  of  the  method  used  in  the  ROSS  model  (Shen  et  al.  1990). 


Lagrangian  discrete-parcel  algorithm 

ta  the  Lagrangian  discrete-parcel  algorithm,  both  the  oil  on  the  water  surface  and  the  oil  in  the  sus- 
pensicsi  are  represented  as  ensembles  of  a  large  number  of  small  parcels.  Each  parcel  has  a  set  of  time- 
dependent  spatial  coordinates  and  a  mass  associated  with  it.  The  movement  of  each  parcel  in  ^'^e  river 
is  affected  by  the  wind,  the  water  current  and  the  concentration  of  surrounding  parcels.  During  each 
time  step,  all  the  oil  parcels  are  first  displaced  according  to  the  current  velocity  and  a  fluctuation  com¬ 
ponent  at  their  respective  locations.  The  turbulent  fluctuation  component,  which  represents  the  effect 
of  horizontal  diffusion,  is  related  to  the  diffusion  coefficient  based  on  the  random-walk  analogy.  After 
all  die  oil  parcels  are  displaced  according  to  advection  and  diffusion,  further  modifications  are  intro¬ 
duced  to  account  for  ail  other  transport  and  weathering  processes.  Due  to  these  processes,  especially 
die  exchanges  between  the  two  layers,  the  mass  of  each  oil  parcel  will  be  different  and  will  change  from 
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time  to  time.  The  number  of  parcels  will  also  increase  with  time  during  the  simulation.  If  a  large  num¬ 
ber  of  parcels  are  released  in  the  river,  and  their  discrete  paths  and  masses  are  followed  and  recorded 
as  functions  of  time  relative  to  a  reference  grid  system  in  space,  then  the  density  distribution  of  th 
particles  in  the  water  body  can  be  interpreted  to  give  the  concentration  of  the  oil. 

The  approach  requires  an  efficient  bookkeeping  procedure  rather  than  the  solution  of  a  large  matrix 
associated  with  a  conventional  Eulerian  finite-difference  or  finite-element  method.  The  algorithm  is 
inherently  stable  with  respect  to  time,  although  the  time  step  should  be  compatible  with  the  grid  size 
and  velocity  for  numerical  accuracy.  Since  the  movement  of  each  parcel  in  the  over  depends  on  the 
distribution  of  the  entire  ensemble  on  the  water  surface  and  in  the  suspension,  all  parcels  must  be 
traced  at  each  time  level  before  proceeding  to  the  next. 

The  detailed  structure  and  implementation  of  the  present  numerical  model  will  be  discussed  and 
presented  later.  In  the  following  sections  the  analytical  formulations  used  for  each  component  of  the 
model  will  be  discussed. 


Model  formulation 


Current  velocity 

Since  the  water  current  affects  advection,  spreading  and  the  exchange  of  oil  between  the  two  layers, 
the  distribution  of  the  magnitude  and  direction  of  the  current  and  the  flow  depth  must  be  determined 
first.  Numerous  numerical  methods  exist  in  the  literature  for  determining  the  two-dimer\sional  flow 
distribution  in  shallow  waters  (Leendretse  1970,  Hamilton  et  al.  1982).  However,  this  type  of  method  is 
time  consuming  and  impractical  for  long  rivers.  A  quasi-two-dimensional  stream-tube  method  is  used 
in  the  present  study.  In  the  present  approach  the  time-dependent  discharge  distribution  Q(x,t)  along 
the  river  is  first  obtciined  through  the  use  of  a  one-dimensional  hydraulic  transient  model,  which  was 
developed  based  on  the  St.  Venant  equatioris: 


dx  dt 


(5) 


dt 


9A 

dx 


+  ^^+gA^-gA(S^-Sf)  =  0 
A  dx  dx 


(6) 


where  x  =  longitudinal  distance  along  the  river 
A  =  flow  cross-sectional  area 
H  =  water  level 
S(,  =  bed  slope 
Sf  =  friction  slope. 

The  ftiction  slope  can  be  calculated  by  Manning's  equation; 
2.21  A^R 


(7) 


in  which  R  is  the  hydraulic  radius  and  n^,  is  the  Manning's  rou^iness  coefficient  of  the  bed.  For  an  ice- 
covered  reach,  the  composite  Manning's  coefficient,  which  accounts  for  the  resistance  of  ice  cover  and 
the  riverbed,  should  be  used  instead  of  nj,.  The  composite  Manning's  coefficient  can  be  Ccdculated  by 
the  Beloktm-Sabaneev  formula  (Larsen  1%9): 


n 


2/3 


(8) 


The  hydraulic  radius  can  be  assumed  to  be  half  of  the  flow  depth  for  an  ice-covered  reach  in  wide 
rivets. 


8 


Once  the  one-dimensional  solution  is  obtained,  the  discharge  Q  can  then  be  distributed  across  the 
width  of  the  river  using  a  stream-tube  model  (Shen  and  Ackermann  1980)  to  give  the  two-dimensional 
velocity  distribution  at  selected  cross  sections. 

For  a  channel  cross  section,  as  shown  in  Figure  4,  the  transverse  distribution  of  the  flow  can  be  deter¬ 
mined  using  the  stream-tube  method.  The  cross  section  is  first  discretized  into  trapezoidal  elements.  By 
applying  Manning's  equation,  the  ratio  of  discharges  between  the  entire  cross  section  and  a  partial 
cross  section  can  be  written  as: 


where  P  =  number  of  trapezoids  in  the  partial  cross  section 

N  =  total  number  of  trapezoids  describing  the  entire  cross  section 
Qp  =  cumulative  discharge  up  to  and  including  the  trapezoid 
Q  =  total  discharge  through  the  entire  cross  section 
Ap,  Rp  =  area  and  hydraulic  radius  of  the  p*  trapezoid,  respectively. 

The  cumulative  discharge  Qp  can  be  computed  by  first  rewriting  eq  9  as: 

p=i 

with 


and 


(9) 


(10) 


(11) 


Qp  -  Qp-i  =  Fq  ApRp^?  (12) 

Based  on  the  computed  distribution  of  Qp,  stream-tube  boundaries  within  the  cross  section  can  be 
determined  by  simple  interpolation.  Once  the  stream-tube  boundaries  are  located,  the  flow  through 
each  stream  tube  is  then  divided  by  the  cross-sectional  area  of  the  stream  tube  to  obtain  the  depth- 
averaged  velocity,  and  the  area  of  the  stream  tube  is  divided  by  its  width  to  obtain  the  width-averaged 
depth.  The  velocity  and  depth  are  then  assigned  to  the  center  of  the  stream  tube.  By  applying  this 
procedure  to  successive  ctoss  sections  along  the  river,  a  two-dimensional  depth-averaged  velocity  and 


a.  Cumulative  discharges  over  a  b.  Area  and  discharge  for  the  Pth  trapezoid, 

partial  cross  section. 

Figun  4.  Method  for  determining  the  transverse  distribution  of  the  flow  in  a  cross  section. 
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Figure  5.  Stream-lube  velocity  distributions  for  a  reach  of  the  Ohio  River. 

width-averaged  depth  distributions  can  be  obtained.  The  direction  of  a  velocity  vector  is  the  same  as 
the  vector  connecting  center  points  of  a  stream  tube  at  successive  cross  sections  in  a  relatively  straight 
reach,  or  perpendicular  to  the  cross  section  in  a  river  bend.  As  an  example,  the  simulated  depth-aver¬ 
aged  velocity  distribution  for  a  reach  of  the  Ohio  River  is  shown  in  Figure  5. 

Once  the  velocity  and  depth  distributions  along  stream  tubes  are  established,  the  distribution  of 
velocity  for  all  the  points  in  a  predefined  grid  system  can  be  obtained  through  linear  interpolations. 
Samples  of  calculated  velocity  and  depth  distributions  are  shown  in  Figures  6  and  7.  The  vedidity  of  the 
stream-tube  method  has  been  examined  by  Shen  and  Ackermann  (1980). 

Adveclion 

During  each  time  step  Af  of  the  oil  spill  simulatim,  the  displacement  AS  of  an  oil  parcel  is  calcu¬ 
lated  by  numerically  integrating  the  drift  velocity  Vt  over  the  time  period  Af.  When  the  Af  value  is 
large,  subintervals  5t^.  are  needed  for  numerical  accuracy.  In  this  case  the  displacement  during  the  time 
interval  Af  is 

AS  =  X^k5fk  (13) 

Jt=l 

where  Vk  =  drift  velocity  of  the  oil  parcel  during  the  time  interval  6fj^ 

AS  =  displacement  during  foe  time  interval  Af 
i  8fk  =  Af. 

The  values  of  6fk  should  satisfy  the  condition  (Roadte  1972,  Qteng  et  al.  1984) 

5fk  S  +  -^1"'  (14) 

Ax  AyJ 

Mfoere  Uk^rid  are  foe  x  and  y  ccHnponents  of  foe  velocity. 

Open  latter  amditim.  In  foe  fnesent  model  foe  advective  velocity  of  each  oil  parcel  is  computed  as: 

V,  «V  (15) 
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where  Vt  is  ttw  drift  velocity  of  an  oil  parcel  and^V  and  V'  are  the  mean  and  turbulent  fluctuation 
compcments  of  the  drift  velocity.  The^mpcaient  V'  is  included  to  simulate  the  horizontal  diffusion  of 
the  oil  parcels.  Hie  formulation  for  V  will  be  discusse<yater.  In  the  suspended  layer  the  mean  drift 
velocity  is  equal  to  the  depth-averaged  current  velocity  Vc,  determined  from  the  stream-tube  simula- 
tkm.  civ  dve  surf»»  layer  ttve  mean  drift  velocity  is  omsidered  to  be  a  vector  surn  of  a  wind-induced 
esad  water-cunent-induced  drift  (Stolzenbadv  et  al.  1977).  The  mean  drift  velocity  Vg  on  the  surface  can 
be  calculated  by  ftte  formula: 

(16) 

»  wind  vdodty  at  10  m  above  ftve  water  surfttce 
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Vc  =  depth-averaged  current  velocity 

=  wind  drift  factor  to  account  for  the  drift  of  the  surface  slick  due  to  wind 
Oj.  =  coiTCCdOTi  factor  for  the  current  on  the  surface. 

A  value  of  0.03  has  been  commonly  used  for  Oy  (Huang  and  Monastero  1982).  The  wind  speed  on  the 
river  is  usually  difterent  ftom  that  at  a  nearby  meteorological  station  (Weeks  and  Dingman  1972).  A 
conecdan  factor  should  be  used  when  the  meteorological  station  is  situated  away  from  the  river. 

Assuming  a  logarithmic  velocity  profile,  the  surface  velocity  of  the  water  current  can  be  related  to 
ttie  depth-averaged  current  velocity  by  the  relationship: 

ttc  =  ^  =  1  +  (17) 

Vc  kVc 

where  Vo  -  surface  velocity 
u»  ss  shearvelodty 
K  s  Kannan  ctmstant  (0.4). 

The  ratio  usually  varies  between  1.1  attd  12. 

hx-caoend  condition.  Oil  spill  movement  tmder  ice  is  a  topic  that  has  received  little  theoretical  or 
laboratmy  tieatxnaat  For  the  suspended  layer  the  trarr^rt  of  oil  droplets  can  be  treated  by  the  advec- 
thre  ditfurion  formulatioiv  as  described  in  eq4.  The  advection  velocity  is  equal  to  the  depth-averaged 
current  vdod^.  The  advection  of  surface  oU  on  die  underside  of  the  ice  cover  is,  however,  more  com¬ 
plicated.  As  suggeded  by  the  experimental  study  of  Cox  and  Schultz  ( 1981),  an  ice  cover  may  be  classi¬ 
fied  into  three  categories  fen'  advection.  The  ice  cover  may  be  smooth,  contain  small  roughness 
dements,  or  contain  korge  roughness  elements.  The  following  discussiext  summarizes  dte  results  of 
their  fimtted  experimental  stiultes  and  serves  as  die  basis  of  the  present  simulatitm. 
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Under  smooth  ice  covers  with  no  current,  the  oil  will  rest  at  an  equilibrium  thickness  that  is  de¬ 
scribed  by  the  empirical  eqxiation 

5«,  =  1.67-8.5(ApJ  (18) 

where  8eq  =  static  equilibrium  slick  thickness  (cm) 

^  =  (Pw  ~  Po)/ Pw'  ttxe  relative  density  difference  between  oil  and  water 
p^y  and  pQ  =  water  and  oil  densities,  respectively  (g/cm^). 

An  ice  cover  is  considered  to  be  smooth  when  the  height  of  the  ice  roughness  is  smaller  than  the  equi¬ 
librium  thickness  of  the  oil,  8gq. 

The  current  velocity  at  which  the  oil  begins  to  move  along  an  ice  cover  is  called  threshold  velocity, 
£ith-  For  a  smooth  cover  the  value  of  was  empirically  determined  to  be  a  function  of  the  oil  viscosity 
Po  and  is  given  as 

Uth  =305.79 /( 88.68 -Po)  (19) 

where  U*  has  the  units  of  cm/s  and  Po  has  units  of  g/cm-s.  Viscosities  for  crude  and  fuel  oils  fall  in  the 
range  of  5-50  centipoises  (1  cp  =  1.0  x  10-2  g/an-s = 2.4  lb/ ft-hr). 

A  rough  ice  cover  has  die  ability  to  retain  the  oil  between  the  roughness  elements.  As  the  current 
velocity  is  increased,  the  oil  will  creep  along  the  upstream  face  of  the  roughness  element  until  it  spills 
over  the  element  and  moves  downstream.  The  threshold  current  velocity  at  which  the  oU  will  move 
downstream  under  a  rough  ice  cover  is  called  the  failure  velocity,  Ufi: 

lifl  =  1.5  { 2  (2^^)  [Oo/^^  (Pw  +  Pb)]  (20) 

where  Oo/w  is  the  oil-water  inteifacial  tension.  The  failure  velocity  Ufi  is  the  velocity  above  which  no  oil 
can  be  contained  upstream  of  a  large  roughness  element. 

If  the  depth-averaged  current  velocity  is  larger  than  the  flireshold  velocity,  iith  or  L/(j,  the  relation¬ 
ship  between  the  current  velocity  and  the  slick  velocity  is  given  as: 

(l-J^f=:  - K -  (21) 

'  ’  0.115  Fi  +  1.105 


with 


F  ^  y  c 

where  V  =  mean  oil  drift  velocity 
Vc  -  current  speed 
K  =  friction  amplification  factor 
Fg  s  slkdcdensimetricFroude  number 
g  -  gravitational  acceleraticm. 


(22) 


The  value  of  K  fa  a  function  of  ttie  roughness  height  of  the  cover.  With  the  limited  data  available,  K  is 
assumed  to  vary  Unholy  between  1.0  for  a  smoofo  cover  and  2.6  for  an  ice  cover  with  Manning's  rou^- 
ness  coefficient  rj  ^  0.0^. 
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Horizontal  diffusion 

The  term  V"  in  eq  15  accovmts  for  the  horizontal  turbulent  diffusion  due  to  the  turbulent  fluctuation 
of  the  drift  velocity.  Based  on  the  random-walk  analysis  (Fischer  et  al.  1979): 

V"  =  (4DT/5f)^/2  (23) 


and 


V'^V'R^e'^'  (24) 

where  5f  =  time  step 

Dj  =  diffusion  coefficient 

Rn  =  normally  distributed  random  niunber  with  a  mean  value  of  0  and  a  standard  deviation 
ofl. 

The  directional  angle  0'  is  assumed  to  be  a  uniformly  distributed  random  angle  ranging  between  0  and  k. 

In  rivers  the  diffusion  coefficient  is  affected  by  the  shear  velocity  u.  and  the  depth  of  flow  h  (Fischer 
etal.  1979): 


Dt=  0.6/1  {J».  (25) 

Hume  experiments  by  Sayre  and  Chang  (1969)  indicated  that  Df  for  surface  dispersants  can  also  be 
calculated  by  eq  25.  Using  Manning's  equation,  we  can  rewrite  eq  25  as: 


DT  =  0.40nbVgl/2;i5/6  (26) 

for  free  surface  flow  (ft^/s)  and 

Dt  =  0.64  (27) 

for  ice<overed  flows  (ft^/s). 

Mechanical  spreading 

The  mechanical  spreading  process  only  affects  the  spreading  of  surface  slicks.  Formulations  used  in 
the  model  will  be  discussed  in  this  section. 

Spreading  in  open  ivater.  A  number  of  theories  have  been  proposed  for  the  process  of  mechanical 
spreading  in  open  waters  (Blokker  1964,  Fay  1969, 1971,  Hoult  and  Suchon  1970,  Mackay  et  al.  1980b). 
In  this  study  Fay's  spreading  theory  (1971)  is  used  because  this  theory  is  based  on  a  rather  comprehen¬ 
sive  descriptian  of  the  spreading  mechanism  and  has  been  verified  by  laboratory  experiments  (Fay 
1971,  Hoult  and  Suchon  1970)  and  other  analytical  solutions  (Farmelop  and  Waldman  1972). 

Fay's  spreading  ttteory  is  derived  for  single-component,  constant-volume  slicks  with  idealized  con¬ 
figurations  in  quiescent  water.  This  theory  considers  the  spreading  of  oil  as  a  result  of  two  driving 
forces — gravity  and  surface  tension — counterbalanced  by  inertia  and  viscous  forces.  The  spreading  of 
an  oil  slick  is  ccsisidered  to  pass  through  three  phases.  In  the  beginning  phase,  orfiy  gravity  and  inertia 
forces  are  important.  In  the  intermediate  phase,  gravity  and  viscous  forces  dominate.  The  firtal  phase  is 
governed  by  tiie  balarKe  between  surface  tension  and  viscous  forces. 

Formulas  for  both  oi\e-dim»isicmal  spreading  and  radial  spreading  at  different  stages  ate  summa¬ 
rized  in  TaUe  1.  The  spreading  rate  during  each  phase  can  be  cbtained  by  taking  time  derivatives  of  the 
formulas  given  in  Table  1 .  The  equations  for  spreading  rates  are  summarized  in  Table  2.  The  time  rate 
of  dhaitge  of  the  oil  volume  in  these  equations  represents  the  change  due  to  weathering  including 
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Table  1.  Spreading  law  for  oil  slicks.  (After  Fay  1971,  Hoult  1979,  Wald- 
man  et  al.  1973.) 


Spreading  Width 

phase _ _ L, 


Gravity-inertia  1.39  (As/lt2)'''3 

Gravity-viscous  1.39  (AjA  2  [3/2  vr-i/2)i/4 

Surface  tension-viscous  1.43  v-*)''''* 

A  =  0.5  volume  of  the  oil  per  unit  length  of  oil  slick. 

A  =  l- 

Pw 

V  =  kinematic  viscosity  of  wafer. 


Radius 


1.14  (AgVt 

0.98  (AxV2  (  3/2  v-l/2)l/b 
1.60(g^t3p;3v-l)l/4 


evaporation  and  vertical  mixing,  and  the  changes  in  oil  volume  distribution  in  various  parts  of  the 
slick.  In  addition  to  the  rates  of  spreading  at  different  phases,  the  times  at  which  phase  transition?  ocdt 
also  need  to  be  determined.  These  transition  times  can  be  obtained  by  letting  equations  in  the  appropri¬ 
ate  phases  equal  each  other  and  solving  for  the  time.  Equations  for  the  transition  times  are  summarized 
in  Table  2. 

Fay  {1969, 1971)  observed  that  the  changes  in  slick  properties  caused  by  weathering  may  result  in 
the  eventual  cessation  of  mechanical  spreading.  Based  on  a  number  of  fields  observations.  Fay  sug¬ 
gested  that 


Af=10^V^^^  (28) 

where  Af  is  the  final  sUck  area  in  m2  and  V  is  the  total  volume  of  the  slick  in  m^.  In  this  study  the  ces¬ 
sation  of  mechanical  spreading  is  considered  to  occur  when  the  slick  thickness  reduces  to  10-5  yw  ni. 

The  formulas  presented  in  Tables  1  and  2  were  derived  for  simple  slick  geometries  that  exist  under 
idealized  conditions.  In  the  simulation  model  the  radial  spreading  formulas  are  used  when  the  slick  is 
nearly  circular,  and  the  one-dimensiortal  formulas  are  used  when  the  slick  area  is  elongated.  A  slick  is 
ccKisidered  to  be  elongated  when  the  aspect  ratio  of  the  slick  area  is  greater  than  3.  The  aspect  ratio 
refers  to  the  length-to-width  ratio  of  the  slick  and  the  orientation  refers  to  the  angle  0  between  the  major 
axis  of  the  slick  and  the  x  axis,  as  shown  in  Figure  8. 

The  orientation  of  fire  slick  is  computed  using  the  moments  of  inertia  and  the  product  of  inertia  of 
the  slick.  The  angle  0  can  be  calculated  using: 

tan  (20)  =  ^^^  (29) 

where 

Z,  =  %2,Iy  =  Zx2P^  =  5:^  (30) 

^4itere  4  and  fy  are  the  moments  of  inertia  of  (he  oil  slick  with  respect  to  the  x  and  y  axes,  respectively, 
and  Pxy  is  the  product  of  inertia.  Once  the  orientation  is  known,  the  x  and  y  coordmates  of  all  the  parti- 
cUs  in  tibe  slide  are  transformed  into  coordinates  in  an  x '  and  y '  coordinate  system  in  which  the  x '  axis 
is  ttte  nve^r  axis  of  the  slick.  The  aqsect  ratio  is  computed  by  the  following  equation,  vdth  x '  and  y ' 
coordinates  of  all  oil  particles: 

Abject  ratio  =  (31) 

2:y' 
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Table  2.  Spread  rates  of  oil  slicks  and  phase  transition  tunes.  (After  Shen 
etal.1990.) 

One-dimensional  slicks,  dLJdt 

Spreading  phase 

Spreading  rate 

Gravity-inertia 

|kio(Ag)i^(L-i/3^f2/3+L2/3r 

■1/3^ 

Gravity-viscous 

it20(Agv-l^r(|U^/8^^,3/8) 

Surface  tension-viscous 

Circular  slicks,  dRJdt 

Spreading  phase 

Spreading  rate 

Gravity-inertia 

Gravity-viscous 

Surface  tension-viscous 

Times  of  phase  transition 

Transition 

One-dimensional  spreading 

Radial  spreading 

Gravity-inertia  to 
gravity-viscous 

(kz.f  /  V 
(k2i)  IvAg/ 

Gravity-viscous  to 
Surface  tension-viscous 

Note: 

=  0.5  volume  of  oil  per  unit  length  along  the  major  axis  of  the  slick  (Lis  a  characteristic 
lengtfi);  k,,)  =  139,  fcjj  =  139,  k^o  =  1.43. 

V  =  Total  volume  of  oil  slick;  fcjj  =  1.14,  k2y  =  0.98,  fcj,  =  1.60. 


pect  ratio  and  orientation. 
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Figure  9.  Division  of  a  nearly  circular  slick  into  pie  - 
shaped  segments. 


Figure  10.  Division  of  an  elongated  slick  into 
strips. 


Oil  Slick  Strips 


This  aspect  ratio  is  used  to  determine  whether  the  spreading  of  the  oil  slick  is  radial  or  one-dimen¬ 
sional.  The  tu>e  of  3  as  the  aspect  ratio  for  transition  from  radial  spreading  to  one-dimensional  spread¬ 
ing  is  subjective  but  gives  reasonable  resrdts.  For  a  nearly  circular  slick  the  slick  area  is  divided  into 
eight  pieshaped  sectors,  as  shown  in  Figure  9.  For  an  elongated  slick  the  entire  slick  is  segmentized 
into  a  series  of  strips,  as  shown  in  Figure  10.  For  both  cases  frie  spreading  rate  of  each  segment  is  calcu¬ 
lated  using  Fay's  formula  shown  in  Table  2,  independent  of  other  segments  in  the  slick.  This  simplifica¬ 
tion  assumes  ttiat  concentration  gradients  at  the  boundaries  between  neighboring  segments  are  negli¬ 
gible  fcff  mechanical  ^reading  (Shen  and  Yapa  1988).  The  segmentizatim  also  allows  for  difr^ent  spread¬ 
ing  rates  in  different  regions  of  the  slick,  thus  providing  a  more  realistic  description  of  the  field  situa- 
tiOTU  For  a  nearly  dicular  slick  the  rate  of  outward  movement  of  an  oil  parcel  along  the  radial  direction 
in  a  particular  pie-shaped  sector  located  at  a  distance  r  from  the  centroid  of  the  slick  can  be  calculated 
from  the  spreading  rate  of  the  mean  radius  7  of  the  sector. 

4L  =  Ldr_.  (32) 

dt  r  dt 

Parcds  scattered  at  distances  far  away  from  the  main  slick  will  be  excluded  from  this  process  since 
small  isolated  patches  of  oil  will  not  be  subjected  to  mechanical  spreading.  In  this  study,  parcels  that 
account  for  the  outer  5%  of  the  total  slick  volume  are  excluded.  This  is  equivalent  to  excluding  the 
parcels  located  at  a  radial  distance  greater  than  2.2  7  from  the  centroid  of  the  slick,  £is  determined  from 
numerical  experiments  for  the  spreading  of  circular  slicks. 

For  an  elongated  slick  tfe  rate  of  outward  movement  in  the  width  direction  of  an  oil  parcel  in  a 
segment  with  mean  lei^fri  I ,  located  at  distance  I  from  the  centroid  of  the  segment,  can  be  calculated 
from  die  rate  of  spreading  of  the  mean  length  of  die  segment  as: 

=  i  (33) 

dt  i  dt  ' 

Parcels  for  die  outer  5%  of  die  dick  vohune  are  not  subjected  to  mechanical  spreading. 

Spreading  under  an  ice  cover.  Hdult  et  al.  (1975)  suggested  that  appreciable  mechanical  spreading  will 
only  oooxr  during  a  oxitbiuous  because  the  oil  thickness  is  stabilized  when  die  equilibrium  thick- 
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ness  for  the  flow  condition  is  reached.  There  is  no  pressure  gradient  or  surface  tension  force  to  cause 
further  spreading  of  the  oil.  The  oil  reaches  an  equilibrium  state  when  cavities  formed  by  the  ice  rough¬ 
ness  contain  a  volume  of  ou  thal  can  decrease  only  with  a  significant  increase  in  the  curren  t  speed.  Since 
a  continuous  spiU  wUl  repeatedly  add  oil  to  jiU  the  cavities,  the  excess  oil  will  effectively  spread  to  tlie 
empty  neighboring  cavities  and  establish  an  equilibrium  state  there. 

The  formula  used  to  model  mechanical  spreading  for  a  continuous  spill  under  ice  is: 

r  =  0.25p.^-^g-  f2/3  (34) 

\  h'  I 

where  r  =  slick  radius 

^  =  (Pw  ~  Po)/Pw-^  the  relative  density  difference 
g  =  gravitational  acceleration 

Qo  =  average  volume  flow  rate  from  the  beginning  of  the  spill 
h  ’  =  half  of  the  root-mean-square  roughness  height  of  the  ice  cover. 

This  equation  is  a  result  of  balancing  the  frictional  drag  from  the  ice  cover  with  the  pressure  drop  that 
occurs  as  the  oil  flows  into  open  cavities.  In  the  simulation  models  it  is  assumed  that  no  mechanical 
spreading  will  occur  for  an  instantaneous  spill  or  after  the  continuous  oil  discharge  stops.  If  the  oil  dis¬ 
charge  is  in  progress  and  the  slick  is  nearly  circular,  the  mechanical  spreading  will  be  calculated  by  eq 
34. 


Shoreline  boundary  conditions 

When  oil  reaches  a  shoreline,  it  will  deposit  along  it.  The  deposited  oil  may  be  re-entrained  into  the 
river  by  the  cument.  The  rate  of  re-entrainment  can  be  formulated  according  to  the  "vulnerability"  in¬ 
dex  of  the  shorelirie  (Gundlach  and  Hayes  1978).  This  index  is  designed  to  reflect  the  enviroirmental 
sensitivity  of  the  shoreline  to  oil  pollution.  For  beaches  of  different  vulnerability  indices,  Torgrimson 
(1984)  suggested  "half-life"  values  to  describe  the  ability  of  the  shore  to  retain  the  oil.  Half-life  is  a  pa¬ 
rameter  that  describes  the  "absorbency"  of  the  shoreline  by  describing  the  rate  of  re-entrainment  of  oil 
after  it  has  landed  at  a  shoreline  location.  Table  3  presents  half-lives  for  different  types  of  shorelines 
along  with  their  vulnerability  indices.  According  to  the  half-life  formulation,  the  fraction  of  the  beached 
oil  re-entrained  into  the  water  body  during  each  time  step  dVj,  is  calculated  from: 


=  1  _  e-kAt  =1-0  5^/^ 

Vb 

where  Vj,  is  the  volume  of  oil  on  the  beach.  The  decay 
ccsnstant  k  can  be  expressed  in  terms  of  the  half-life  X  as: 


(35) 

Table  3.  Shoreline  descriptor  and  default 
parameters. 


k  =  [-In  (1/2)]  A.  (36) 

For  the  purpose  of  calculatir^  evaporation  rate,  the 
beached  oil  is  assumed  to  occupy  a  predefined  width  of 
each  computaticHial  grid  along  the  ^oreline.  This  width 
is  determined  by  the  user  according  to  the  field  corwii- 
tiom. 

Evaporation  and  dissolution 

Evf^xmtion.  Evaporation  is  die  main  process  diat  caus¬ 
es  die  kMS  in  cril  volume  durir^  weathering.  Mackay  et 
al.  (l^K)a)  developed  a  fcmnuk  for  calculating  the  rate 


Shoreline  descriptor 

Half-life 

(hr) 

Vulnerability 

index 

Imtantaneous  rejection* 

0 

0 

Exposed  headUmd 

1 

1 

Wave-cut  platform 

1 

2 

Pocket  beach 

24 

3 

Sand  beach 

24 

4 

Sand  and  gravel  beach 

24 

5 

Sand  and  cobble  beach 

8760 

6 

Exposed  tide  fiats 

1 

7 

Sheltered  rock  shore 

8760 

8 

Sheltered  tide  flats 

8760 

9 

Sheltered  marsh 

8760 

10 

Land 

8760 

11 

*  This  option  b  added  in  the  present  computer  model 
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of  evaporation.  The  volume  fraction  of  oil  evaporated  is  determined  as 


F  =aE(l/C)[lnPo  +  ln(CKEf +  I/P0)] 


(37) 


where 


E  =  K^t 
Ke 

Km 

A 

V 

R 

T 

v; 


modification  coefficient  introduced  in  the  present  study  to  account  for  the  effect  of  the 
change  in  oil  properties  on  the  evaporation  rate  due  to  emulsification  and  other  weather¬ 
ing  processes 

evaporative  exposure  term,  which  varies  with  time  and  environmental  conditions 
Kf^AvARTV'j 

mass  transfer  coefficient  (m  /s)  (0.0025  U  w^a) 

spill  area  (m^) 

molar  volume  (m^/mole) 

gas  constant  [82  x  10-6  m^/fmole-K)] 

surface  temperature  of  the  oil  (K),  which  is  generally  close  to  the  ambient  air  tempera¬ 
ture  Tg 
volume  (m3). 


The  initial  vapor  pressure  Pg  in  atm  at  the  temperature  Te  is  determined  from 


lnPo=10.6[l-To/TE]  (38) 

where  Tq  is  the  mitial  boiling  point  (K).  The  constant  C  can  be  determined  by  the  relationship  Tp_C  = 
constant.  Values  of  C  for  =  283  K  are  given  in  Table  4.  For  crude  oils  of  different  API  index  values,  C 
at  Te  =  283  K  are  given  in  Table  5  along  with  Tg.  This  table  can  be  replaced  by  the  following  functional 
relationships  obtained  through  curve  fitting  (Shen  et  al.  1990); 

C  =1158.9 API (39) 

and 

To  =  542.6  -  30.275API  -t-  l.SfiSAPl^  +  34.39AP/3  -i-  0.0002604  AP/^  (40) 

The  API  index  and  the  specific  gravity  of  the  oil  are  related  by 

Specific  gravity  =  141 5/  (API  131 5) .  (41 ) 


Table  4  Suggested  evaporation  parameters  for 
various  petroleum  fractions  (Tg  =  283  K). 
(After  Mackay  et  al.  1980b.) 


Table  5.  Suggested  evaporation 
parameters  for  various  crude  oils 
(Te  =  283  K).  (After  Mackay  et  ai. 
1980.) 


Oil  type 

TJK) 

C 

P„  (atm) 

Specific  gravity 
(gcitr^) 

Motor  gasoline 

Summer 

314 

5.99 

0313 

API 

C 

P. 

Winter 

308 

6.23 

0.39 

10 

1.0 

89.2 

366 

0.044 

Aviation  gasoline 

341 

2.81 

0.12 

12 

0.986 

69.4 

348 

0.088 

Diesel  fuel 

496 

5S7 

3.4x10-* 

15 

0.966 

52.1 

339 

0.13 

Jet  fuel 

418 

5.06 

6.0  X 10-® 

20 

0.934 

34.7 

329 

0.18 

No.  2  funutce  oil 

465 

7.88 

1.1  X 10^ 

25 

0.904 

27.2 

330 

0.17 

Lube  (heavy  and  light) 

583 

8.61 

132  X  1(H 

30 

0.876 

22.33 

325 

0.21 

ffeavy  gascril 

633 

8.99 

2.0  xKH 

35 

0.850 

19.5 

314 

031 

Residuals 

783 

337 

735  X  10-» 

40 

0.825 

17.9 

304 

0.45 

L^htgascui 

473 

637 

8.1  X 10-* 

45 

0.802 

16.4 

283 

1.004 

19 
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Monongahela-AUegheny  river  system. 


The  molar  volume  of  oil  is  required  in  eq  37.  The  value  ‘A'  uie  molar  volume  can  vary  bet\vt*en  1 50  x 
KH  and  600  x  10^  m^/mole.  Typical  values  of  molecular  weight  of  various  oil  components  are  given 
by  Moore  et  al.  (1973).  The  compositions  of  several  different  oils  are  given  along  with  the  concentration 
and  molecular  weight  Mj  of  each  component.  Based  on  these  the  molecular  volume  can  be  computed 
(Shen  et  al.  1990).  The  volume  Vq  for  each  time  step  is  determined  by 

Vo'  ^  Vo  +  AV„  (42) 

where  V^  =  effective  initial  volume  of  surface  oil 

Vq  =  iiutial  volume  of  oil  spilled  on  the  water  surface 
AVn,  =  volume  increase  in  the  surface  layer  due  to  the  exchange  between  the  two  layers. 

Dissolution.  Cohen  et  al.  (1980)  developed  a  method  to  calculate  the  total  dissolution  rate  N: 

N^KA^S  (43) 

where  N  =  total  dissolution  rate  of  the  slick  (g/hr) 

K  =  dissolution  mass  transfer  coefficient,  assumed  to  be  0.01  m /hr 
As  =  slick  area  (m2) 

S  =  oil  solubility  in  water. 

Huang  and  Monastero  (1982)  suggested  that  for  a  typical  oil  the  solubility'  can  be  calculated  by 


S  = 


(44) 


where  is  the  solubility  for  fresh  oil  and  t  is  time  in  hours.  Huang  and  Monastero  suggested  a  typical 
value  of  30  g/m^  for  S^.  In  the  present  study  the  formula  of  Lu  and  Polak  (1973)  is  used,  which  pro¬ 
vides  more  infonnation  on  solubility.  In  this  method  the  rate  of  dissolution  is  calculated  as 

Td  (45) 


where  is  the  rate  of  dissolution  (mg/m2-day).  For  three  oil  samples  tested  the  coefficient  c  and  d  are 

given  in  Table  6. 


Table  6.  Oissolution  coefficients  at  25‘'C. 
(After  Lu  and  Polak  1973.) 


OiV 

type 

API 

c  d 

(mg/m^)  d/day} 

KS,-cd 
{Vg^  hr) 

No.  2  fuel  oil 

35.5 

1043 

0.423 

0.0184 

Crude  oil 

38.6 

8915 

2.380 

0.884 

Bunker  C  oil 

14.8 

459 

0.503 

0.0104 

Mixing,  emulsification  and  sedimentation 

Tenns  representing  mass  exchanges  in  the  surface  and  tfie  suspended  layers  in  eq  1  and  4  are  calcu¬ 
lated  at  the  end  of  each  time  step.  The  voliune  of  oil  mixing  into  the  suspended  layer  per  time  step  in  a 
grid  is  calculated  as  "yCgAxAyAf.  Similarly  the  volume  of  suspended  oil  resurfaced  is  calculated  as 
<sCyV\^AyM.  The  volume  of  oil  settled  to  the  channel  bed  is  calculated  as  pCViArAyAf . 

COMPUTER  MODEL  AND  CASE  STUDY 

Based  on  tihe  analytical  formulation  presented  in  the  previous  chapter,  we  developed  a  computer 
model  for  tt»<Dhio-MQnc«E^3hela--Allegheny  River  (OMA)  system  (Fig.  11).  This  chapter  describes  the 
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no. 

Origiu;-  of  looil 
coordinatei  ,y,,> 

i/t’\’rc'db 

River 

1 

0.0, 0.0 

0.0 

Ohio 

2 

348200.0,  1388IX).0 

-90.0 

Ohio 

3 

378400.0,  -42200.0 

-90.0 

Ohio 

4 

392000.0,  -84800.0 

-180.0 

Mononghel.! 

5 

315800.0,  -104800.0 

-180.0 

Monunghcla 

6 

383200.0,  -36800.0 

360.0 

Allcghuny 

7 

592000.0,  -768000.0 

90.0 

Alli'ghviiy 

model  along  with  a  case  study.  It  should  be  not-  Table?.  Relationships  between  local  and  global 

ed  that  the  following  discussions  deal  with  coordinate  systems, 
many  programming  details.  For  a  user  to  ap¬ 
ply  the  model  to  the  OMA  system,  a  menu-driv¬ 
en  model  (Yapa  et  al.  1990a)  provides  a  simple 
means  for  spill  simulations. 

Model  implementation 

The  grid  system 

In  the  computer  simulation  a  two-dimen¬ 
sional  reference  grid  is  used  to  discretize  the  riv¬ 
er.  This  reference  grid  is  required  to  identify  the 
locations  and  velocities  of  oil  parcels.  The  grid 
boxes  in  this  model  are  squares  of  equal  size. 

The  location  of  each  grid  box  is  identified  by  its 
X  and  y  ndices. 

The  reference  grid  described  is  adequate  for 
rivers  with  only  one  orientation,  e.g.  from  east 
to  west,  or  with  only  small  variations  in  their 
orientation.  Unfortunately  the  OMA  river  sys¬ 
tem  has  a  complicated  configuration  with 
changing  orientations.  Seven  sets  of  local  coor¬ 
dinates  cire  introduced,  as  shown  in  Figure  12, 
for  different  regions  in  the  system  to  accommo¬ 
date  the  changing  orientations.  During  the  sim¬ 
ulation,  transfonnation  between  a  local  coordi¬ 
nate  system  ana  the  global  system  is  required. 

The  local  x-axis  is  set  to  be  in  the  main  flow  di¬ 
rection  toward  upstream.  Relatior\ships  be¬ 
tween  the  global  and  local  coordinate  systems 
are  given  in  Table  7. 

Besides  the  fixed-grid  systems  a  temporary 
moving  grid  system  traveling  with  die  slick  is 
used  in  the  model  to  save  memory  space.  This 
moving  grid  system  will  just  be  long  enough  to 
cover  the  length  of  the  oil  slick  and  the  travel 
distance  in  one  time  step  of  the  Bow  computa¬ 
tion.  In  a  long  simulation  the  slick  may  become 
too  long  to  be  covered  by  the  moving  grid  because  of  bank  deposition  and  rejection.  For  this  case  there 
are  two  cations,  one  of  which  is  to  make  the  length  of  the  moving  grid  long  enough  to  cover  the  longest 
slick  during  die  spill  or  the  total  length  of  the  river.  However,  this  option  is  limited  by  the  memory 
space  of  microcomputers.  Hence,  the  second  option  is  used.  In  this  option  the  model  selects  a  reason¬ 
able  length  of  moving  grid,  both  for  computer  memory  and  for  the  length  of  the  oil  slick,  for  example, 
30  miles  long.  The  moving  grid  will  follow  the  slick's  leading  edge  and  cover  a  30-mile-long  segment  of 
the  river  upstream  from  the  leading  edge.  The  oil  outside  the  moving  grid,  if  it  exists,  will  generally  be 


I  and  X  Ohio  River 
X.ICandX  Monongahela  River 
3land3ac  Allegheny  River 


Figure  12.  Local  and  global  coordinate  sj,  Jems. 


A  serks  of  successive  river  bends  can  cause  difficulty  in  defining  the  river  boundary.  In  Figure  13 
there  are  four  continuous  river  bends.  There  are  more  than  three  sets  of  boundary  points  that  intersect 
wiffi  a  kxal  /  grid  liiw.  hi  ffiis  model  all  land  boundary  points,  except  for  the  two  extreme  end  points. 
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Local  X  or  i 


Figure  13.  Method  for  determining  the  land  grid  along  a  series  of  river  bends. 

are  treated  as  island  botindary  points.  The  boundary  grids  along  a  series  of  bends  can  be  divided  into 
three  types,  as  shown  in  Rgure  13.  They  arc  island  grids,  bank  grids  and  their  combinations. 

Release  of  oil 

The  oil  released  into  the  river  is  represented  by  a  user-specified  number  of  particles.  Each  particle 
can  be  considered  as  a  parcel  of  oil.  The  user  can  select  a  particle  number  up  to  1000  for  each  layer  at  the 
beginning  of  the  simulation.  During  the  simulation  the  number  may  increase  to  3000  or  even  more  than 
5000  due  to  the  exchanges  between  the  surface  and  the  suspended  layers.  For  simulation  accuracy, 
approximately  500  particles  in  each  layer  is  suggested  to  be  used  at  the  beginning  of  the  simulation. 
However,  in  an  emergency  situation  the  user  may  reduce  the  particle  number  to  50  or  less  to  reduce  the 
computing  time.  The  volume  of  oil  released  into  each  layer  is  determined  by  the  user's  input.  Parcels  in 
the  same  layer  have  the  same  volume  at  the  beginning.  The  total  number  of  parcels  in  each  layer  is  the 
same  at  the  beginning,  so  that  the  oil  is  released  at  a  constant  rate  in  each  layer.  For  example,  consider  a 
spill  of  10,000  gallons,  with  7000  gallons  in  the  surface  layer  and  3000  gallons  in  the  suspended  layer.  If 
the  user  chooses  1000  parcels  in  each  layer,  then  each  surface  oil  parcel  will  contain  7  gallons,  and  a 
suspended  one  will  contain  3  gallons  at  the  beginning. 

An  oil  spin  will  be  treated  either  as  an  instantaneous  spill  or  a  continuous  spill.  When  the  spUl  dura¬ 
tion  is  zero,  the  spill  will  be  treated  as  instantaneous.  Otherwise,  it  will  be  treated  as  a  continuous  spill. 

Velocity  and  depth  distributions 

Analjrsis  of  the  transport  of  oil  spill  parcels  in  a  river  requires  well-defined  water  velocity  and  depth 
distributions.  To  calculate  tttese  distributions,  various  information  is  needed,  including  ihe  total  dis¬ 
charge,  river  stage,  cross-section  geometry,  shoreline  characteristics,  ice  condition  and  wind  condition. 
The  velocity  and  depdt  distributions  in  the  river  are  calculated  in  two  stages.  First,  a  one-dimensional 
unready  flow  model  is  used  to  generate  dischar^  and  river  stages  at  specified  cross  sections  referred 
to  as  no^.  Rtr  the  CMA  model  the  Flowsed  model  is  used  for  this  purpose  {Johnson  1982).  The  oil 
^dil  model  uses  this  information  to  compute  two-dimensional  velocity  and  depth  distributions. 

In  the  oil  ^>ill  model,  branches  are  used  to  describe  the  channel  configuration.  The  beginning  and 
end  of  these  brandies  are  generally  selected  to  coindde  wifii  nodes  of  the  one-dimensional  unsteady 
flow  model.  In  cases  where  it  is  rwcessary  to  have  more  branches  in  the  oil  spill  model  than  in  the 
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a.  Nodes.  b.  Branches. 


Figure  14.  Schemalization  of  a  river  reach  with  an  island. 

vmsteady  flow  model,  water  levels  at  intermediate  points  are  obtained  through  interpolation.  When  an 
island  is  encoimtered,  there  will  be  branches  on  each  side  of  the  island  (Fig.  14).  Each  of  these  branches 
includes  a  portion  of  the  stream  tubes.  The  information  available  is  the  magnitude  and  direction,  i.e. 
the  X  and  y  components,  of  velocities  at  stream-tube  centers  at  all  cross  sections. 

The  assignment  of  velocity  and  depth  through  the  entire  river  requires  the  establishment  of  a  grid 
system,  as  shown  in  Figure  15.  Velocities  and  depths  obtained  from  the  stream-tube  analysis  are  inter¬ 
polated  to  all  grid  boxes.  River  boundaries  for  this  grid  system  are  defined  by  boundary  grid  boxes.  For 
each  grid  along  the  x-direction,  two  corresponding  y-grid  boxes  are  used  to  define  the  location  of  the 
riverbanks.  In  places  where  there  are  islands,  more  y-grid  boxes  are  used  to  define  the  upper  and  lower 
boundaries  of  islands.  The  model  can  handle  any  number  of  islands  in  the  river.  However,  if  more 
islands  are  intersected  by  a  vertical  line  along  which  the  x  coordinate  is  constant,  more  space  will  be 
needed  for  the  input  data,  and  simulation  speed  will  be  slowed  down.  Defirung  islands  in  the  grid 
system  is  independent  of  defining  islands  in  a  branch  configuration.  This  allows  for  a  greater  flexibility 
in  the  simulation  of  oil  slick  transformation.  For  example,  a  small  island  that  covers  orUy  three  or  four 
grid  boxes,  which  may  be  too  small  to  be  defined  in  the  branch  corifiguration,  can  be  easily  defined  in 
the  grid  box  system.  All  islands  must  be  defined  in  the  grid  box  system. 


1  2  3  4  5  6  7  8  9  tO  11  1213  14  15  16  17  18  19  20  2122  23  24  25 

i  •  grids 

Figure  15.  Grid  boxes  and  river  boundary  representation. 
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The  velocity  and  depth  within  a  box  are  assumed  to  be  uniform.  Only  boxes  found  within  the  de¬ 
fined  river  boundaries  are  assigned  velocities  and  depths.  Generally  spacing  between  cross  sections  in 
the  stream-tube  computation  is  much  larger  than  the  grid  size  Ax.  The  following  procedure  is  adopted 
for  assigning  velocity  and  depth  to  grid  boxes  based  on  the  velocity  and  depth  at  cross  sections.  In  the 
following  discussion  the  term  "velocity/ depth  point"  is  used  to  represent  a  point  at  which  a  computed 
velocity/ depth  is  located.  At  first  cdl  velocity  and  depth  points  computed  from  the  stream-tube  analy¬ 
sis  will  be  assigned  to  the  boxes  within  which  they  lie.  If  more  than  one  velocity  and  depth  point  fail 
into  a  box,  the  average  value  of  velocity/ depth  is  assigned  to  that  box.  In  the  next  step,  more  velocity / 
depth  points  will  be  obtained  through  interpolation  between  two  successive  cross  sections.  The  num¬ 
ber  of  interpolation  points  between  two  successive  cross  sections  can  be  changed  through  user-defined 
input.  These  velocity/ depth  points  obtained  through  interpolation  will  now  be  assigned  to  grid  boxes 
using  the  same  method  as  before.  After  this,  if  there  is  any  box  without  an  assigned  velocity /depth,  an 
average  velocity/ depth  calculated  based  on  neighboring  boxes  will  be  assigned.  Upon  the  completion 
of  these  steps,  the  entire  river  is  scaimed  to  ensure  that  velocity/depth  points  are  assigned  to  all  boxes. 

Oil  slick  orientation 

The  surface  slick  can  spread  either  axisymmetrically  or  one-dimensionally.  The  determination  of 
which  type  of  spreading  prevails  depends  on  the  shape  of  the  slick,  i.e.  the  aspect  ratio.  Therefore,  after 
the  oil  particles  moved  to  their  new  positions,  the  centroid,  aspect  ratio  and  orientation  of  the  oil  slick 
are  calculated.  These  calculations  will  consider  only  oil  particles  in  the  water  grids. 


Table  8.  Bound¬ 
ary  types  of  locks 
and  dams. 


Code 

no. 


Rejection 
_ rate—. 


13 

14 

15 

16 
-IL. 


Figure  16.  Index  system  far  locks  and  dams. 


Lock  and  Dam 

1 3  »  Particle  moves  opposite  local  x-coordmate. 

14  =  Partide  moves  along  local  >^axis 

15  »  Particle  moves  opposite  local  y-axis. 

16  «  Particle  moves  along  local  x-axis. 


Boundary  conditions  and  locks  and  dams 

Twelve  half-life  values  were  defined  in  Table  3  for  shoreline  boundaries.  Oil  slicks  passing  a  lock 
and  dam  will  be  completely  mixed  into  die  water  column.  The  program  is  directed  to  this  condition 
thrmi^  the  use  of  ^oreline  boundary  codes.  Four  new  boundary  condition  codes  are  added.  As  de¬ 
fined  in  Table  8,  code  numbos  13-16  refer  to  lock  and  dams  with  four  different  orientations  (Fig.  16) .  A 
set  of  special  rejection  rate  codes  are  given  for  diese  conditions  to  direct  the  program  to  weU-mixed 
conditions.  Oil  parcels  deposited  <m  the  banks  at  a  lock  and  dam  are  instantly  rejected  into  the  water, 
lids  is  adiieved  by  using  a  set  of  shoreline  codes  0  and  6  as  in  Table  8.  This  set  of  codes  can  also  be  used 
at  ottier  dhorelir^  locations  for  instantaneous  rejection  if  needed. 
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START 


Figure  17.  The  computer  model  ROSS2. 


Model  structure 

Hie  comfniter  code  is  written  in  standard  Fortran  77.  No  compiler-spedfic  extensions  have  been 
used.  Hie  code  has  been  tested  only  on  32>bit  machines  with  Fortran  77  compilers  that  support  32-bit 
cOQ^tati<xis.  Hie  program  can  be  run  on  other  machines  (e.g.  24-bit)  with  minor  modifications.  Ap¬ 
pendix  B  gives  a  Usting  of  all  the  subrcmtines  used.  Hw  model  is  named  ROSS2,  and  the  model  struc¬ 
ture  is  shown  in  Hgure  17.  AH  variables  are  dimensioned  to  have  sufficient  storage  to  run  for  each 
s^pent  of  die  river  systexn.  Hie  velocity  and  dqpth  distributicHis  wlQ  be  kept  only  for  a  reach  long 
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enough  for  the  oil  slick  to  travel  one  un-  Table  9.  Cross  references  of  subroutines, 

steady  flow  calculation  time  step.  The 
model  can  be  run  for  a  new  spill  simula¬ 
tion  or  continue  from  a  previously  ter¬ 
minated  run  by  changing  weather  and 
oil  characteristics. 

Subroutine  cross  reference 

The  model  consists  of  the  main  pro¬ 
gram  ROSS2  and  37  subroutines,  as  sum¬ 
marized  in  Table  9.  An  alphabetically  ar¬ 
ranged  listing  of  the  subroutines  is  giv¬ 
en  in  Appendix  B. 

The  common  block  includes  arrays 
and  their  definitions  shared  by  subrou¬ 
tines.  Sizes  or  arrays  can  be  adjusted  in 
this  block  according  to  the  moving  grid 
size,  maximum  particle  number,  river 
segment  number,  maximum  island 
number  and  stream-tube  number.  Table 
10  gives  definitions  of  important  vari¬ 
ables  included  in  the  common  block. 

When  the  model  is  applied  to  other  riv¬ 
ers,  dimensions  should  be  adjusted.  Oth¬ 
er  variables  used  in  the  program  are  pre¬ 
sented  in  Appendix  A. 

Main  program 

The  main  program  is  used  to  estab¬ 
lish  the  link  between  subroutines,  be¬ 
sides  interactions  going  through  the 
common  block.  The  major  functions  of 
ttie  main  program  can  be  suirunarized 
as  follows: 

•Read  oil  spill  information  in  the  data 
filezzxsp/. 

•Read  boundary  grid  box  numbers 
and  set  boimdary  conditions  from 
xxx.igr  and  xxx.bnd. 

•Call  con  to  read  back  the  simulation  results  if  it  continues  a  job  from  a  previously  terminated  run. 

•CaU  ckout  to  adjust  the  moving  grid  position  based  on  the  oil  slick  position. 

•Call  veldep  to  calculate  velocity  and  depth  distributions  of  the  river  reach  superimposed  by  moving 
grids  during  every  unsteady  flow  time  step. 

•Call  presfis  and  presur  to  release  oil  particles  at  equal  intervals  for  the  oil  spill  duration. 

•CaU  advsus  and  advsur  to  advect  oU  particles  in  both  the  suspended  layer  and  the  surface  layer. 

•CaU  emdam  to  mix  oU  particles  into  the  water  column  when  passing  locks  and  dams. 

•CaU  orient  to  calculate  the  oU  stick  orientation,  determining  whether  the  slick  is  quasi-one-dimen- 
sional  or  quasi-axisymmetric;  caU  pievol  to  calculate  the  volume  of  each  pie  segment  for  an  axisym- 
metric  stick;  or  <aU  strvol  to  calculate  the  volume  in  each  smaU  strip  for  a  one-dimensional  slick.  AU 
of  Uiese  are  for  medhanical  spreading  calculations. 

•CaU  emulfii  and  redepn  to  cxrmpute  oU  volumes  into  the  water  column,  resurfaced  to  the  water  sur- 


Program 

Calk 

Input  files 

ross2 

input,  ckout,  positn,  veldep,  ohndc, 
ohvad,  advsur,  presus,  advsur,  advsus, 
bndloc,  checkl,  ckout,  con,  rejioc, 
emdam,  orient,  pievol,  strvol,  emulfn, 
redepn,  aread,  evapor,  disolu,  lumper, 
exclutn,  bndrej,  shift,  sprdax,  sprdld, 
rearrg,  ploter 

River  division  data, 
spill  data, 

shoreline  condition, 
ice  regions 

advsur 

gauss,  randu,  bndrej,  rejioc 

None 

advsus 

gauss,  randu,  bndrej,  rejioc 

None 

aread 

None 

.None 

bndloc 

None 

None 

bndrej 

rejioc 

None 

cakij 

None 

None 

caldm 

None 

None 

dieckl 

cakij 

None 

ckout 

caldjxaklm 

None 

con 

conchk 

Previous  simulation 
results 

conchk 

cakij 

Previous  simulation 
results 

disolu 

None 

None 

emdam 

None 

None 

emulfn 

None 

None 

evapor 

None 

None 

exchan 

None 

None 

input 

None 

Spill  constants, 
river  divison  data 

lumper 

None 

None 

lumperl 

None 

None 

ohndc 

None 

Flow  data 

ohvad 

None 

None 

pievol 

None 

None 

ploter 

None 

None 

posit! 

None 

None 

presur 

advsur,  gauss,  randu,  bndrej,  rejioc 

None 

presus 

advsus,  gauss,  randu,  bndrej,  rejioc 

None 

rearrg 

lumperl 

None 

relca 

None 

None 

reloctt 

None 

None 

reloxy 

None 

None 

relxy 

None 

None 

shift 

None 

None 

sprdax 

None 

None 

sprdld 

None 

None 

strvol 

None 

None 

veldep 

ohndc,  ohvad 

Geometry  data, 
flow  data 
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Table  10.  Definition  of  variables  in  the  common  block  (dimension  parameters). 


Variable  Value _ Definition 


numbm 

37 

Maximum  number  of  branches  of  river  segments 

numdep 

66 

Maximum  number  of  sounding  depth  locations  at  all  cross  sections 

numice 

10 

Maximum  number  of  ice  regions 

numpar 

15000 

Maximum  number  of  oil  particles 

numpie 

8 

Number  of  pies  into  whirii  the  axisymmetrical  slick  is  divided 

numvlO 

2 

Number  of  positions  where  the  velocities  are  assigned  zero 

nuinpos 

26000 

Maximum  number  of  moving  grids 

numsec 

330 

Maximum  number  of  cross  sections  among  river  segments 

numstb 

11 

Number  of  stream  tubes 

nwth 

10 

Number  of  shorelines,  including  islands,  locks  and  dams 

nudxx 

1100 

Maximum  grid  number  in  x  direction  of  moving  grids 

numxgd 

4600 

Total  number  of  grids  in  the  x  direction 

nseg 

5 

Maximum  number  of  river  segments 

nisi 

4 

Maximum  number  of  islands  a  cross  section  can  meet 

nout 

800 

Maximum  number  of  particles  out  of  moving  grids  that  can  be  plotted 

face  and  deposited  onto  the  riverbed  during  the  time  step,  as  well  as  the  cumulative  values. 

•Call  aread  to  calculate  the  slick  area  in  water,  under  ice  and  on  land  for  the  calculation  of  evapora¬ 
tion  and  dissolution. 

•Call  evapor  and  disolu  to  calculate  evaporation  and  dissolution. 

•Call  lumper  and  exchan  to  calculate  oil  exchanges  between  the  suspended  layer  and  the  surface 
layer. 

•CaH  sprdax  for  an  axisymmetrical  oil  sUck,  or  call  sprdld  for  a  one-dimensional  slick,  to  calculate 
mechanical  spreading. 

•Call  bndrej  to  reject  oil  particles  back  into  the  river  according  to  the  bank  retarding  condition  (half- 
life). 

•Call  rearrg  to  adjust  particle  sizes,  Le.  lump  small  particles  together  and  divide  large  particles  to 
satisfy  the  requirements  of  the  random-walk  simulation  and  the  computer  memory  space. 

•Call  plater  to  calculate  the  one-dimensional  cross-sectional  average  oil  concentration  distribution 
along  the  river. 

•Save  the  simulation  results  into  corresponding  output  files  at  a  desigjnated  output  frequency. 
Subroutines 

Since  aU  subroutines  and  related  variables  are  listed  in  Appendices  A  and  B,  and  some  simple  sub¬ 
routines  can  be  easily  understood  with  comments  in  the  listing,  no  further  explanation  will  be  given. 
Block  diagrams  for  some  main  subroutines  are  given  in  Figure  18  to  assist  in  understanding  the  pro¬ 
gram. 

Case  stilly 

The  model  is  applied  to  the  oil  spill  in  the  Monongahela  River  on  2  January  1988  (Miklaucic  and 
Saseen  1989).  Simulatkm  results  for  tiie  leading  edge  and  concentration  distribution  of  tire  oil  plume 
will  be  compared  with  the  field  data. 

Model  parameters 

A  number  of  empirical  parameters,  which  can  affect  the  simulation  results,  are  to  be  selected  by  the 
user.  In  tiius  section  fire  effects  of  these  parameters  will  be  discussed  to  assist  tire  user  in  making  appro¬ 
priate  selections. 

Qrte  (rf  file  parameters  Lhat  will  affect  the  movement  of  the  slick  is  the  shoreline  boundary  ccmditicm. 
Iftider  an  instantaneous-rejection  shoreline  boundary  condition,  oil  will  not  deposit  on  riverbanks.  Par- 
tides  moved  to  the  shme  will  be  reacted  to  water  instantly.  This  wiU  give  a  shorter  slick,  and  its  leading 
edge  will  advance  at  a  faster  rate. 

nwoe^iatify  fitc  simuiation  should  not  be  affected  by  the  spill  time  step.  However,  numerically  a 


28 


Figure  18.  SuimniHms. 


END 


START 
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e.  Particle  size  adjustment. 
Figure  18  (cant’d). 


large  time  step  can  lead  to  inaccurate  results,  especially  for  the  non-instantaneous-rqection  condition. 
In  the  model,  boundary  rejection,  exchange  between  layers,  mechanical  spreading,  evaporation  and 
dissolution  are  done  once  every  spill  time  step  to  save  computation  time  and  memory.  The  spill  time 
step  can  affect  the  slick  length  through  all  of  these  simulation  components.  If  the  time  step  is  too  large, 
discontinuous  distributions  of  tire  portions  of  oil  deposited  on  shoreline,  on  the  water  surface  and  in 
the  suspoision  will  appear.  When  the  instantaneous-rejection  condition  is  used,  the  effect  of  the  spill 
time  step  is  small 

The  random-walk  theory  requires  a  large  number  of  particles  with  equal  mass.  Numerical  experi¬ 
ments  show  that  the  equal  mass  requirement  can  be  relaxed  if  the  particle  number  is  large  enough.  In 
ttie  present  model  the  particle  munber  increases  when  the  particle  mass  becomes  more  tmeven  due  to 
the  0cchange  between  the  surhice  and  suspended  layers  The  model  also  keep>s  at  least  one  particle  in 
eadi  grid  that  contains  oil 

To  demcmstrate  tfie  effects  of  the  ^oreline  conditirm,  simulations  for  the  time  step  At  and  the  initial 
paitide  number  are  carried  out  for  the  1988  oil  spill  in  the  Monongahela  River  for  a  duration  of  10 
days.  The  model  parameters  used  in  these  simulations  are: 

•  Bur^ant  velocity;  r\, = 0.01  ft/s; 

•  Resurfacing  oodffikient:  a  =  1.0; 

•  Emulsification  coeffidenb  ^  =  lO-^/s; 

•  Hiveibed  deposition  coeffidenfc  10-^/s; 
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•  Initial  particle  niimben  N^,  =  40,  or  400  each  layer; 

•  Spill  time  step:  At  =  5, 15  or  30  minutes;  and 

•  Evaporation  modification  coefficient;  %  =  03. 

Ice  covers  can  affect  the  movement  of  oil  in  the  river.  According  to  the  model  the  advection  velocity 
of  the  oil  slick  under  a  smooth  ice  cover  will  only  be  reduced  by  10-15%.  The  ice  cover  has  very  little 
effect  in  the  advection  of  suspended  oil.  For  the  1988  spill  the  air  temperature  was  below  freezing,  as 
shown  in  Figure  19.  Ice  may  have  existed  in  the  Ohio  River  during  the  later  part  of  the  simulation 
period.  However,  since  the  ice  covers  formed  in  the  Ohio  River  are  generally  smooth,  the  existence  of 
an  ice  cover  will  be  neglected  in  the  sample  simulations  for  simplicity. 


Figure  19.  Air  temperature  variation  during  the  simulation  period, 
Pittsburgh,  Pennsylvania. 


Figure  20  shows  simulation  results  with  the  idealized  instantaneous-rejection  shoreline  condition 
and  a  partial-rejection  shoreline  condition  determined  for  the  OMA  system.  Table  11  and  Figures  21 
and  22  compare  simulated  and  observed  locations  of  fire  leading  edge  of  the  slick  and  the  concentra¬ 
tion  peedc.  These  results  show  fiiat  the  slick  locations  simulated  using  the  instantaneous-rejection  shore¬ 
line  condition  do  not  agree  with  the  observed  locations,  alttiough  the  errors  are  relatively  small.  Simu¬ 
lations  using  the  partial-rejection  shoreline  condition  compared  well  with  the  field  observation.  Simu- 


Table  11.  Observed  and  simulated  locations  of  the  leading  edge  and  the  concentration  peak.* 


Leading  edge  Concentration  peak 

(rtoermile)  (rwer  mile) 

Date  Observed  Simulated  Observed  Simulated 


Instantaneous  shoreline  rejection 


2  Jan 

-25.0 

-25.0 

-25.0 

-25.0 

3  Jan 

0.0 

-1.21 

-43 

-234 

4  Jan 

25.0 

22A4 

no 

19.04 

5  Jan 

45.0 

5^14 

no 

49.90 

6  Jan 

65+ 

7432 

no 

7238 

7Sr 

7432 

60.0 

7238 

7Jan 

73* 

9839 

65.0 

9637 

88** 

9639 

77.0 

%.07 

8Jan 

89.0+ 

1153 

773 

112.9 

99J0r 

1153 

84.0 

11Z9 

9  Jan 

940* 

125.1 

no 

1223 

108" 

125.1 

86.0 

1223 

10  Jan 

112.0 

1253+ 

96.0 

125.1+ 

Leading  edge 
(rwer  mile) 

(river  mile) 

Date 

2  Jan 

Partial  shoreline  rejection 
-25.0  -25.0  -25.0 

-25.0 

3  Jan 

0.0 

0.41 

-4.5 

-Zll 

4  Jan 

25.0 

21.06 

no 

14.06 

5  Jan 

45.0 

43.11 

no 

3234 

6  Jan 

65* 

58.82 

60.0 

52.43 

78" 

58.82 

60.0 

52.43 

7  Jan 

73* 

77.62 

65.0 

64.06 

88" 

77.62 

77.0 

64.06 

8  Jan 

89.0* 

90.43 

77.0 

78.41 

99.0" 

90.43 

84.0 

78.41 

9  Jan 

94.0* 

1013 

86.0 

8539 

108" 

1013 

86.0 

85.39 

10  Jan 

112.0 

109.7 

96.0 

9331 

*  Initial  petttide  nundier  in  each  layer  is  400;  yill  time  step  is  900  s. 

*  Vieoty  and  AMcs-Kedzkna  (1989). 

**Beilceyel8La989). 
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Figure  20.  Simulation  results 
shamng  the  variation  of  con¬ 
centration  ^MrSmthn. 
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Figure  21.  Comparison  of  simulated  and  observed  leading  edge  positions;  =  400,  At  =  300  s. 


latkxns  uang  small  particle  numbers  generally  resulted  in  shorter  slicks  with  higher  concentration.  The 
efiect  of  time  step  At  is  shown  in  Kgure  23.  This  figure  shows  that  the  slick  slowed  down  with  the 
increase  in  At.  However,  this  etiect  is  small. 

A  summary  of  computing  time  used  for  all  of  the  sample  simulations  run  on  a  SUN  SPAC  stationl+ 
is  givm  in  Table  12.  This  t  •  showc;  that  the  computing  time  decreases  with  a  reduction  in  particle 
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At,  Spill  Tima  Step  (min.) 

Figure  23.  Effect  of  time  step  At  on  leading  edge  positions. 

Table  12.  Comparison  of  computing  time. 

Parameters _ CPU  time  (s)  Parameters _ CPU  time  (s) 

Instantaneous  rejection  i-utial  rejection 

No  =  400;  A/ =  300  s  23,855  N,  =  400;A.  =300  s  67A54 

No  =  400;  At  =  900  s  10,797  N^  =  400;  A/  =  900  s  32,248 

No  =  400;  At  =  1800  s  3,612  Np  =  400;  At  =  1800  s  19,258 

number  and  an  increase  in  id.  The  use  of  the  instantaneoits-rejection  shoreline  condition  reduces  the 
computing  time  by  approximately  60%.  It  is  recommended  that,  for  field  applications  in  theOMA  river 
system,  five  partial-rejection  shoreline  cwidition  with  At  =  900  s  and  =  400  be  used,  considering  the 
simulation  accuracy  and  computing  time.  The  ratio  between  the  real  time  and  the  computing  time  is 
about  1:12.  This  ratio  is  smaller  at  the  early  stages  of  simulation.  For  emergency  applications  an  instan- 
taneous-rejectian  shoreline  condition  along  with  a  smaller  N„  and  a  larger  At  ^ould  be  used.  The  ratio 
between  the  real  time  and  the  computing  time  is  about  1:240  when  Af  =  30  minutes  and  =  40  are 
used. 

Applications  to  non-oil  spills 

As  mentioned  earlier  the  model  can  be  applied  to  non-oil  spills.  Non-oil  spills  can  be  simulated  by 
setting  appropriate  values  for  foe  spill  parameters  described  earlier.  For  neutrally  buoyant  miscible 
materials,  the  buoyant  velocity  should  be  set  to  zero. 

INPUTDATAFILES 

There  are  two  categodes  of  inputs  that  are  required  to  run  the  model.  The  first  category  can  be  con- 
Sklered  to  be  fixed  data  for  a  given  river  reach.  This  is  foe  information  required  to  describe  foe  shore¬ 
line  ami  crow-sectional  geometry  of  foe  river.  Normally  there  is  no  need  to  adjust  foese  data.  The  sec¬ 
ond  category  Includes  various  parameters  foat  may  be  adjusted  according  to  the  flow  condition  and 
dharaeteitetks.  To  completdy  understand  foe  setup  of  foe  data  files,  it  is  helpful  to  go  forough  foe 
8ie|>4y-stiqf>fmtoeciuie  foatfolloi^  Sample  input  data  files  are  given  later.  If  the  user  is  only  interested 
to  ac^ueth^  pararneteis,  ciumging  the  spill  kxation  or  establishing  riew  stage  and  clischarge  ccmditicms 
OR  a  iMevknisly  litodded  river,  no  new  data  files  need  be  created. 
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Data  file  creation 

For  demonstration  purposes  the  Ohio-Monongahela  river  combination  will  be  used  as  an  example. 
This  river  shows  the  complexity  of  a  river  system  that  Uie  model  can  handle.  Where  the  data  files  for 
the  Ohio-Monongahela  river  combination  are  not  sufficient  to  describe  the  detail,  other  river  data  will 
be  presented.  Eight  data  files  exist  for  inputting  information  into  the  computer  model; 
xxx.geo  River  geometry,  cross  sections  cind  branch  geometry 
xxxAgr  Boundary  grid  box  numbers 
xxx.ice  Ice  regions,  ice  roughness  and  oil  viscosity 

xxx.flo  Water  level  and  discharge  at  nodes  determined  from  the  unsteady  flow  model 
xxx.bnd  Index  for  half-life  value  assignments  to  shore  grids 
xxx.spl  Oil  parameters,  spill  location,  wind  condition  and  ambient  temperatirre. 
xxx.ndc  Information  to  transform  computed  unsteady  flow  results  to  what  are  needed  in  the  oil  spill 
model  flow  computation 

xxx.msc  River  segment  definitions,  parameters  for  exchange  between  two  layers  and  information  need¬ 
ed  to  coordimte  transfonnation  data. 

The  first  step  in  preparing  the  data  file  is  to  make  a  sketch  similar  to  Figure  24  to  describe  the  branch 
and  cross-section  numbering  system  used  in  the  data  file.  The  procedure  for  obtaining  this  sketch  is  as 
follows: 

1.  Etetermine  the  number  of  river  segments.  In  this  model  the  entire  Ohio-Monongahela-AUegheny 
river  system  is  divided  into  seven  segments,  as  described  earlier.  The  numbering  of  segments 
depends  on  the  river  combination.  The  Ohio-Monongahela  has  five  river  segments  numbered 
from  downstream  to  upstream.  Each  segment  has  its  own  local  coordinate  system,  which  can  be 
transformed  to  the  global  coordinate  system  using  the  transformation  data  in  data  file  ohmon.int. 
All  geometry  data  will  be  prepared  segment  by  segment  with  independent  numbering.  To  ensure 
that  aU  water  grids  are  covered,  cross-section  numbers  will  have  overlaps  in  the  vidnity  of  the 
division  line  between  river  segments. 

2.  Etetermine  the  number  of  branches.  Each  branch  must  contain  at  least  two  unsteady  flow  nodes 
and  one  stream-tube  cross  section. 

3.  Number  the  cross  sections  in  consecutive  order  from  the  downstream  end  to  the  upstream  end  of 
the  segment  Around  islands,  first  number  the  bottom  side  up  to  the  last  cross  section  prior  to  the 
confluence  into  a  single  channel,  then  go  back  to  the  top  side  and  continue  the  numbering  se¬ 
quence. 


Figure  24.  D^nition  sketch  of  the  river  numbering  system. 


37 


4.  Detennine  the  number  of  stream  tubes  in  each  branch.  The  total  number  of  stream  tubes  always 
remains  the  same.  The  stream  tubes  must  be  divided  according  to  the  ratio  of  the  flow  split  around 
the  island. 

Once  the  branches  and  cross  sections  for  the  river  are  established,  scaled  maps  of  the  river  and  cross 
sections  are  used  to: 

1.  Establish  the  global  x-y  Cartesian  coordinates,  to  be  superimposed  over  the  river. 

2.  Establish  the  reference  coordinates  on  the  left  bank  (looking  downstream). 

3.  Ehgitize  the  entire  river  boundary  and  cross-section  locations  simultaneously. 

4.  Digitize  the  cross-section  geometries  by  measuring  a  transverse  distance  from  the  reference  point 
on  the  left  bank  to  a  corresponding  sounding  point. 

5.  Determine  the  way  to  assign  flow  direction  to  each  cross  section.  For  the  cross  sections  along  river 
bends,  use  1  as  the  direction  index,  which  implies  that  the  flow  directions  will  be  perpendicular  to 
the  cross  section.  Along  straight  reaches,  especially  in  reaches  with  islands,  0  is  chosen  as  the  flow 
direction  index.  This  implies  that  the  flow  direction  will  be  along  the  line  connecting  correspond¬ 
ing  stream-tube  centers  in  two  neighboring  cross  sections. 

6.  Schematize  boundary  boxes,  i.e.  for  every  x-grid  there  exists  corresponding  upper  and  lower  ri- 
verbanks  and  island  boundary  boxes. 

The  files  are  generally  broken  up  into  blocks  and  cards.  A  block  covers  a  group  of  data  that  may 
ccmtain  one  or  more  card  types.  A  card  type  is  one  line  of  specific  data,  which  is  sometimes  repeated. 
(For  example.  Block  4  inxxx.geo  has  Cards  1  and  2,  where  Card  2  is  repeated  as  many  times  as  needed.) 
By  inspecting  file  example  of  a  card  and  comparing  it  to  the  complete  sample  data  set,  it  is  easy  to  see 
how  the  entire  file  comes  toother. 

Most  of  the  data  read  into  the  model  are  in  list-directed  I/O  (free  format).  If  column  numbers  are 
shown,  foe  data  must  be  formatted  accordingly.  Otherwise  it  is  necessary  to  put  orUy  one  space  or 
comma  between  each  number  in  a  card. 

In  foe  following  sectiorrs,  data  files  of  foe  Cfoio-MOTiongahela  river  combination  will  be  used  as 
examples  to  explain  foe  setup  of  data  cards. 

Geometry  data:  ohmortgeo 

Ohmon.gm  consists  of  five  segments  of  geometry  data  for  foe  five  corresponding  river  segments. 
Each  segment  has  four  blocks.  None  of  foe  data  in  this  file  needs  to  be  adjusted  from  one  spill  to  foe 
r«xt.  The  user  may  choose  to  add  additional  cross  sections,  change  foe  number  of  branches,  etc.  Howev¬ 
er,  the  user  is  cautkmed  to  consult  this  section  b^e  making  such  changes.  Since  all  segments  have  foe  same 
format,  only  one  will  be  explained  here. 

ohmon.gar.  Hock  1,  branch  and  grid  mfomuition 

gffldl 

Example: 

24  16 

VariaHe  Typeand  OAumn 

name _ bmgth _ number _ Description _ 

nbtndi  Integer  -  Number  of  branches 

kintm  btteger  -  Number  of  velocity  interpolations  between  cross  sections 

in  a  stream  tube 
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Card  2  (one  number  for  each  branch) 

Example: 

12  22  32  44  45  64  89  104  105  121  131  138  139  141  150  208  210  218  228  233 

Variable  Type  and  Column 

name _ length _ number _ Description _ 

lcstsq(i)  Integer  -  Last  cross  section  in  each  branch.  In  the  last  branch,  use  the 

second  to  last  cross  section.  There  must  be  a  total  of  nbranch 
data.  If  the  line  is  not  long  enough,  continue  to  the  next  line 

ohmon.geo:  Block  2,  cross  section  and  connection  information 

Card  1  (one  card  for  each  cross  section) 

Example; 

1  16840.66  7741.73  -0.44720  0.89443  11  11  2  0  0.0  0 


Variable 

name 

Type  and 
length 

Column 

number 

Description 

j 

Integer 

— 

Cross-section  number 

x 

Real 

— 

X  coordinate  (in  global  coordinate  system)  of  reference  point 
on  left  bank 

y 

Real 

— 

y  coordinate  (in  global  coordinate  system)  of  reference  point 
on  left  bank 

cst 

Real 

— * 

Cosine  value  that  the  angle  cross  section  makes  with  posi¬ 
tive  x-axis  (global) 

snt 

Real 

Sine  value  that  the  angle  cross  section  makes  with  positive 
x-axis  (global) 

nstube(i) 

Integer 

- 

Numb«  of  stream  tubes  at  current  cross  section 

numcon(i) 

Integer 

If  aU  stream  tubes  continue  to  next  cross  section  undivided, 
=  11;  if  stream  tubes  divided  into  two  channels  from  a  main 
charmel,  =  12;  if  stream  tubes  from  this  channel  and  another 
channel  connect  to  next  section  foat  is  in  foe  main  channel,  = 
21 

Next  cross  section  connecting  to  foe  current  cross  section. 
For  a  divided  charmel  around  an  island,  this  represents  foe 
first  cross  secticm  connected  to  in  foe  lower  division  from 
foe  main  channel  cross  section 

r\fuco(i) 

Integer 

nseco{i) 

Integer 

For  a  divided  channel  around  an  island,  this  represents  the 
first  cross  section  connected  to  foe  upper  division  from  foe 
main  channel  cross  section  (if  no  isl^d,  =  0).  If  this  is  foe 
first  cross  section  in  the  upper  branch,  =  888.  If  this  is  foe  last 
in  the  upper  branch,  =  999. 

stub(i) 

Real 

* 

If  a  bramh  has  small  islands,  =  discharge  fraction  cm  left  side 
of  the  island;  if  no  island,  =  0.0 

ibend(i) 

Integer 

— 

Flow  directicm  index;  for  cross  sections  in  a  river  bend,  =  1; 
otherwise  =  0 

ohmon.geo:  Block  3,  cross-section  geometry 
Cardl 


Example: 
6  39 


Variable 

mme 

Type  and 
length 

Column 

number 

Description 

j 

Integer 

— 

Cross-section  number  (for  checking) 

nslsctQ) 

Integer 

Number  of  sounding  depths  (points  along  the  bed  in  the 
transverse  direction)  used  to  describe  the  cross-section  ge¬ 
ometry 

Card  2  fas  manv  cards  as  required  to  input  all  sets  of  Ywid-zl 

Example: 

0.0  660.0 

540.0  660.0 

540.0 

5%.0  ... 

Variable 

Type  and 

Column 

mme 

length 

number 

Description 

ywid(i,j) 

Real 

- 

Transverse  distance  from  the  reference  shore  to  the  sound¬ 

ing  depth  in  i*  cross  section 

2(y) 

Real 

Bed  elevation  at  the  transverse  point  (sounding  depth) 
for  the)**'  transverse  point  (sounding  depth)  for  the ;'**'  cross 

_ section _ 

Note:  Block  3  must  be  repeated  lcstsq(nbmch)  times  (Le.  number  of  cross  sections  defined). 

ohmon.geo:  Block  4,  defines  all  grid  boxes  that  are  to  be  assigned  zero  velocity 


Cardl 

Example: 

0 


Variable 

name 

Type  and 
length 

Column 

number 

Description 

nzrvb 

Integer 

- 

Number  of  boxes  to  assign  zero  velocities 

^2 

VariaUe 

name 

Type  and 
Imgth 

Column 

number 

Description 

izibx 

Integer 

Integer 

mm 

x-grid  number  of  i*  box  to  have  zero  velocity 
y-grid  number  of  i*  box  to  have  zero  velocity 

Note:  There  are  must  be  imvb  pairs  of  izrbx(i)  and  iziby(i).  Data  can  be  continued  to  as  many  lines 
asneeded. 


For  file  od^  river  segmmt,  just  repeat  block  1  to  block  4.  This  is  the  end  of  data  file  ohnwn.geo. 
ke  candid  data:  ckmoxdee 

TbechmmJoe  file  contains infonnation  identifying  ice  regions.  An  ice  regicsi  is  defined  by  a  range  of 
grid  ham  fiutf  are  ioe-ooveted;  tsing  (he  example  in  Figure  25  an  ice  region  may  be  identified  as  ex¬ 
tending  fixim  grid  (14^)  to  grid  (18,7).  The  ke  regkm  then  covets  every  grid  from  (143)  to  (18,7)  as 
diomi  in  (he  figure  This  infonnation  is  used  when  determmmg  if  spreading  arid  advection  takes  place 
under  an  ke  cover  or  on  open  water. 


ohmonAce:  Block  1 

Cardl 

Example: 

0035  125 


Variable 

name 

Type  and 
length 

Column 

number 

Description 

anice 

atniuo 

Real 

Real 

- 

Manning's  n  for  ice  roughness 

Viscosity  of  oil  (gm-cm/s)  (poise) 

Card  2 
Example: 

2 

Variable 

name 

Type  and 
length 

Column 

number 

Description 

nkei^ 

Integer 

- 

Total  number  of  ice  regions 

Example: 
14  8  18 

7 

Variable 

name 

Type  and 
ImgOt 

Column 

number 

Description 

nioexl(i) 

nioeyl(i) 

nioex2(i) 

t<key2ii) 

integer 

Integer 

biteger 

biteger 

— 

X  grid  at  file  beginning  of  ice  region 
y  grid  at  the  begiiuiing  of  ice  region 

X  grid  at  the  end  of  ice  region 
y  grid  at  the  end  of  ice  region 

Hus  is  ttie  end  of  data  file  oftmon.ics. 


chinonJo 

Die  chmmflo  file  contains  the  unsteidy  flow  time  step,  the  water  level  and  the  discharge  at  ^ch 
node  in  file  river  as  defined  by  the  one-diznemional  unst^y  flow  model  (Yang  et  al.  1990).  Also  in- 
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eluded  are  the  ice  conditions  for  each  cross  section  in  the  river.  These  data  are  separated  from  the  ice 
region  data  in  ohmonJee.  The  oil  spUl  simulation  model  converts  this  information  into  boundary  condi¬ 
tions  for  each  river  branch. 

This  file  consists  of  three  blocks  of  information.  All  blocks  are  listed  below  with  descriptions  and 
corresponding  components.  Blocks  2  and  3  must  be  repeated  every  time  when  the  velocities  and  depths 
are  updated  in  the  model,  i.e.  every  time  step  of  the  one-dimerisional  flow  model.  Therefore,  the  data  in 
this  file  need  to  be  adjusted  on  a  more  regular  basis. 

If  the  format  or  column  numbers  are  shown,  the  data  must  be  formatted  accordingly.  Otherwise  it 
necessary  to  have  only  one  space  or  comma  between  the  data. 

For  the  Ohio-Monongahela-AUegheny  river  system,  the  difference  in  branch  configuration  between 
the  oO  spill  model  and  the  unsteady  flow  model  is  built  into  the  oil  spUl  model  through  subroutine 
ohndc.  The  number  of  branches  in  this  data  file  is  the  same  as  the  number  of  branches  in  the  unsteady 
flow  model.  The  discharges  and  water  levels  at  the  nodes  of  the  branches  of  this  data  file  should  be 
interpolated  to  the  cross  sections  of  the  oil  spill  model.  If  the  two  models  have  the  same  number  of 
branches,  tiien  subroutine  ohndc  should  be  changed  to  the  following  three  lines,  which  are  already  in 
the  program.  These  lines  are  presently  given  as  comments,  which  can  easily  be  activated: 
do230i=l,nsect 
read(7,*)  wlsct(i),  qsct(i) 

230 continue 

ohmonflo:  Block  1,  time  step  in  one-dimensional  model 

C-acd-i 

Example: 

6.0 

Variable  Type  and  Column 

name _ length _ number _ Description  _ _ 

ufdt  Real  -  Unsteady  flow  time  step  (hr) 

Card  2  (Comments  line) 

Example 

branch  qu  qd  hu  hd  region  1 1/1/88  hour=0.0 

These  are  the  variable  name,  date  and  time  of  flie  data  to  be  given  in  the  following  block. 

Blodc  2,  discharge  and  water  level 

l(one  card  for  (me  branch) 

Example: 

1  52810.  52810,  625.19  6251)0 

2  51761.  52810.  625.25  625.19 


Variabk 

nami 

Typeand 

lengOt 

Cdumn 

number 

Description 

1 

boega 

— 

Branch  number  for  checking 

qu(0 

Real 

- 

Discharge  at  upstream  0ul  of  each  branch  (c£s) 

Real 

- 

Discharge  at  downstream  end  of  each  aiKh  (cfo) 

WIJ(0 

Beal 

- 

Water  level  at  upstream  end  of  each  branch  (ft) 

wd(i) 

Real 

“ 

Water  level  at  upstream  end  of  each  branch  (ft) 

Hbir.  Card  1  needs  to  be  sweated  ib  many  times  as  tite  number  of  brandies  of  this  river  segment. 


42 


ohmon.flo:  Block  3,  ice  conditions  (thickness  information) 

Card  1 
Example: 

1 

Variable 
name 

idnfo 


Card  2 
Example: 

2  open 

Variable  Type  and  Column 

name  length _ number _ Description _ 

is  i4  -  ID  number  of  the  cross  section  with  ice-covered  conditions 

word  a4  -  Cross-section  ice  cover  condition,  "full"  =  fully  covered, 

"part"  =  partially  covered,  "open"  =  open  water.  If  word  = 
"full,"  then  the  card  has  only  one  value  and  that  is  the  ice 
thickness  across  the  river  for  that  x-section.  If  word  =  "part," 
then  card  3  must  have  an  ice  thickness  defined  at  each  verti¬ 
cal  line  defining  the  x-section  =  nslsct(is}.  If  all  numbers  don't 
_ fit  into  one  card,  as  many  cards  as  necessary  may  be  u*^. 

Card  2  (for  fully  covered  cross  section) 

Card  3  must  foUow  every  card  2  if  word  =  "full"  or  "part."  If  word  =  "open,"  card  3  is  not  needed. 

Example: 

0.6 

Variable  Type  and  Column 

name  Imgth  number  Description 


fulti 

Real 

— 

Ice  thickness  (ft)  of  hilly  covered  cross  section.  Only  one  val¬ 
ue  is  read  as  input,  and  it  will  be  assigned  to  the  entire  cross 
section 

Card  3  (for  partially  covered  cross  section) 

Example: 

03  0.4  05  05  0.4  0.2 

Variable 

name 

Type  and 
Imgth 

Column 

number 

Description 

tice(i,}) 

Real 

- 

Ice  thickness  (ft)  of  partially  covered  cross  section.  There  must 
be  one  value  for  each  sounding  depth  of  die  cross  section 

with  one  additional  value  for  the  extreme  left  of  the  cross 
secticNn  where  depth  sounding  is  not  input  through  data  be- 
_ cause  it  is  always  assumed  to  be  zero _ 

Note:  Block  2  followed  by  3  must  be  repeated  for  every  one-dimensional  unsteady  flow  time  step. 
This  is  the  csid  of  data  file  ohmon.^. 


Type  and  Column 

length  number  Description 

Integer  -  Number  of  cross  sections  with  ice-cover»>d  condition.  If  no 

ice-covered  section  exist,  set  icinfo  =1  and  then  in  card  2,  de¬ 
fine  the  section  to  be  open 


Boundary  grid  data:  ohmonigr 

okmonAgr.  Block  1,  boundary  grid  boxes 

Card  1 
Example: 

343  121  130  125  125 


Variable 

name 

Type  and 
length 

Column 

number 

Description 

1 

i4 

x-grid  box  number 

igrilb(i) 

15 

— 

y-directiem  grid  box  number  of  lower  river  boundary  for  i* 
x-grid  (water-side  grid) 

igriub(i) 

i5 

y-direction  grid  box  number  of  upper  river  boundary  for  fl' 
x-grid  (water-side  grid) 

igrlb(j/i) 

i5 

y-direction  grid  box  number  of  lower  boundary  of  island 

for  i*  x-grid  (land-side  grid) 

igrub(j,i) 

i5 

“ 

y-direction  grid  box  number  of  upper  boundary  of  island 
for  i*  x-grid  (land-side  grid) 

Note:  Block  1  must  repeat  ntogridx  (no.  of  grids  in  x-direction)  times. 


This  is  the  end  of  data  file  ohmonAgr. 

Shoreline  boundary  condition  data:  ohmon.bnd 

The  dhmon.bnd  file  consists  of  three  blocks  of  information  All  blocks  are  listed  below  with  compo¬ 
nents  and  descriptiort  Most  of  the  data  read  into  the  model  are  in  free  format.  If  the  format  or  column 
nuntbcrs  are  shown,  the  data  must  be  formatted  accordingly.  Otherwise  it  is  necessary  to  have  only 
one  space  or  comma  between  the  data.  All  files  with  the  bnd  extension  follow  the  same  format. 

ohmon.bnd:  Block  1,  half-life  index. 

Card  2  (oiie  card  for  each  range  of  grid  boxes) 

Example: 

1  1  81  1 


Variable 

name 

Type  and 
length 

Column 

number 

Description 

k 

Integer 

- 

Store  number.  Lowery  (river)  =  1;  upper y  (river)  =  2;  lower 
y  (first  island)  =  3;  upper  y  (first  island)  =  4;  etc. 

lfix»n 

Integer 

Beginning  limit  (grid  box  number)  for  half-life  designation 
to  shore 

Ito 

Integer 

— 

Ending  limit  (grid  box  number)  for  half-life  des^nation  to 
^ore 

icode 

btteger 

btteger  identifying  which  of  the  16  half-life  values  are  to  be 
assigned  to  a  grid;  for  determinirtg  the  value  for  lock  and 
dams,  refer  to  Figures  15  and  16 

Card  last  (must  be  included) 

0  0  0  0 

This  is  file  end  erf  data  file 

OOi^datafdkmcKispi 

the  file  conHSts  two  blocks  of  infonnatkm  with  a  vatyii^  or  constant  nvunber  of  cards 

ineach  blodc.  AD  blocks  am  fisted  bdow  wifii  the  components  and  description.  Most  of  the  data  read 
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into  the  model  are  in  free  format.  If  the  format  or  column  numbers  are  shown,  the  data  must  be  format¬ 
ted  accordingly.  Otherwise  it  is  necessary  to  have  only  one  space  or  comma  between  the  data.  All  files 
with  the  spl  extension  follow  the  same  format. 

ohmon.sph  Block  1,  Oil  spill  and  simulation  parameters. 

Card  1  (river  name — for  identification  only) 

Example: 

Ohio  Monongahela 


Variable 

name 

Type  and 
length 

Column 

number 

Description 

text 

a20 

- 

Any  text  to  identify  the  river  system 

Card  2  (Type  of  oil.  for  identification  only) 

Example: 

Fuel  No.  2 

Variable 

Type  and 

Column 

name 

length 

number 

Description 

fueltp 

a20 

- 

Text  for  identifying  the  oU  type 

C?ir(d3 

Example: 

24.0  4  1 

1  1  0 

1800.0  - 

-1.0 

Variable 

Type  and 

Column 

name 

leng& 

number 

Description 

totime 

Real 

_ 

Time  of  oil  spill  simulation  (hr) 

ievery 

foteger 

Frequency  of  obtairring  output  from  simulation  results  from 
ohmon  and  plotter.  For  example,  if  ievery  =  4,  output  is  writ¬ 
ten  every  four  time  steps 

ipotl 

Integer 

Two  optitMis:  one  (1)  results  in  output  of  fixed  data,  such  as 
cross^ecticm  geometry  data  and  ^oreline  conditions;  zero 
(0)  results  in  no  output 

ipot2 

Integer 

Two  opticms:  one  (1)  results  in  output  of  computed  veloci¬ 
ties  to  a  data  file  for  plotting;  zero  (0)  results  in  no  printout 

ipot3 

Integer 

Two  options;  one  (1)  results  in  output  of  particle  location  to 
a  data  file  for  plotting;  zero  (0)  results  in  no  printout 

ipot4 

Integer 

Two  options:  one  (1)  results  in  output  of  number  plot  of  par¬ 
ticle  distribution  to  a  data  file  for  plotting;  zero  (0)  results  in 
no  printout 

spltim 

Real 

Duration  of  oil  spill  in  seconds.  If  =  0,  the  spill  will  be  treat¬ 
ed  as  instantaneous;  otherwise  it  will  be  continuous 

difiur 

Real 

Horizontal  diffusirsi  coefficient  (ftV  s)  for  river.  If  the  defoult 
formulation  is  derired,  set  this  value  to  -1.0 
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Card  4 
Example: 

400  705000.  900.  .84  1.4116-5  7.55e^  1.14  1.1  1.6  1.4  1.4  1.4 


Variable 

Type  and 

Column 

name 

length 

number 

Description 

ntotal 

Integer 

- 

Number  of  parcels  defined  in  each  layer  at  the  beginning  of  the  sim¬ 
ulation  (currently  400) 

spvol 

Real 

- 

Total  volume  of  oil  spilled  (gal.) 

spgoii 

Real 

- 

Specific  gravity  of  oil 

aniu 

Real 

- 

I^ematic  viscosity  of  water  (ft^/ s) 

sigma 

Real 

- 

Surface  tension  of  oil  (lb/ ft) 

ak2i 

Real 

- 

Fay's  gravity-inertia  phase  spreading  coefficient  (axisymmetrical) 

al<2v 

Real 

- 

Fay's  gravity-viscous  phase  spreading  coefficient  (axisymmetrical) 

ak2t 

Real 

Fay's  surface  tension-viscous  spreading  coefficient  (axisymmetri¬ 
cal) 

akclO 

Real 

Fay's  or  Waldman's  gravity-inertia  spreading  phase  coefficient  (one¬ 
dimensional) 

akc20 

Real 

- 

Gravity-viscous  phase  spreading  coefficient  (one-dimensional) 

akcSO 

Real 

Surface  tension-viscous  phase  spreading  coefficient  (one-dimension¬ 
al) 

Card  5 
Example; 

342650. 

-97800.  .7063e-2  .1873e-27.88  465.0 

Variable 

Type  and 

Column 

name 

Imgth 

number 

Description 

spx 

Real 

- 

x<oordinate  of  spiU  site  (ft) 

spy 

Real 

- 

y-coordinate  of  spill  site  (ft) 

vmuni 

Real 

- 

Molar  volume  of  oil  (ft^/ mol) 

soluni 

Real 

- 

Solubility  of  fresh  oil  (lb/ ft^) 

cevp 

Real 

- 

Coefficient  C  of  evaporation  characteristics  of  oil 

toevp 

Real 

Boiling  point  temperature  of  oil  (K).  If  you  define  a  value  of  less 
than  1.0  for  toevp,  the  program  will  define  the  evaporation  charac¬ 
teristics,  using  fitted  curves.  Therefore,  the  input  values  of  cevp  and 
toevp  have  no  influence  on  computations,  although  they  are  read 

okmonspl:  Block  2,  wind  and  temperature  data. 

Card  1  (1  card  for  each  time  step 

Example: 
02  11.80 

2150  2400 

Variabk 

Type  and 

Column 

name 

Imgth 

number 

Description 

vwc 

Real 

- 

Empirical  coefficient  for  modifying  the  wirul  speed  on  land  to  that 
on  river 

vwmag 

Real 

- 

Wind  speed  (ft/s) 

-  «  - 

tilBlSI 

Real 

•m. 

Wind  direction  (angle  measured  clockwise  from  the  north  in  de¬ 
grees) 

t«vf 

Real 

Air  temperature  (®F) 

Mole:  C!^  5  must  be  repeated  for  every  tune  step  of  dte  simulation. 


IMs  is  tfte  end  d  data  file 
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Bridging  data  betioeen  the  unsteady  flow  model  and  the  oil  spill  model:  ohinon.ndc 

The  ohmon.ndc  file  consists  of  two  blocks  of  information  with  a  varying  or  nonvarying  number  of 
cards  in  each  block.  All  blocks  are  listed  below  with  the  components  and  description.  Most  of  the  data 
read  into  the  model  are  in  free  format.  If  tine  format  or  column  numbers  are  shown,  the  data  must  be 
formatted  accordingly.  Otherwise  it  is  necessary  to  have  only  one  space  or  comma  between  the  data. 
AH  files  with  the  ndc  extension  follow  the  same  format. 

ohmon.ndc:  Block  1,  stream-tube  cross-section  distribution 

Card  1  (Comments) 

Example: 

Region  1:  part  1 


Variable 

name 

Type  and 
length 

Column 

number 

Description 

dummy 

a4 

Dummy  variable  for  skipping  a  line  in  the  data  file 

Card  2 
Example: 

5  1  45 

Variable 

name 

Type  and 
length 

Column 

number 

Description 

j 

Integer 

- 

Branch  number  in  this  river  segment 

ns 

Integer 

- 

Number  of  stream-tube  cross  sections  in  diis  branch 

jar(j) 

Integer 

- 

List  of  numbers  of  all  stream-tube  cross  sections  in  the  branch 

Note:  Card  2  must  repeat  for  aU  branches  in  this  segment. 


ohman-ndc:  Block  2,  location  of  stream-tube  cross  sections 

Example: 

2  2.11  1  0.09 


Variable 

name 

Type  and 
Imgth 

Column 

number 

Description 

j 

Integer 

_ 

Cross-section  number 

ihup 

Integer 

- 

Number  of  upstream  unsteady  flow  node 

xup 

Real 

- 

Distance  between  cross  section  and  upstream  node 

ihdp 

Integer 

- 

Number  of  downstream  unsteady  flow  node 

xdp 

Real 

- 

Distance  between  cross  section  and  downstream  node 

Note:  2  must  repeat  as  many  times  as  the  number  of  cross  sections  in  this  river  segmmt. 


Repeat  blocks  1  and  2  as  many  times  as  the  number  of  river  segments.  This  is  the  end  of  data  hie 
obmon.ndc. 


Miscethmeous  data:  ahtrununsc 

The  chmorumsc  61e  consists  of  two  blocks  of  information  widi  a  varying  or  nonvarying  nurnb^  of 
cards  in  eadt  Idock.  All  blocks  ate  listed  below  widi  the  components  and  desaiption.  Most  of  die  data 
read  into  the  model  ate  in  bee  format  If  the  format  or  column  numbers  are  shown,  the  data  must  be 
ficamatted  aoocnding^y.  Odterwise  it  is  necessary  to  have  only  one  sp>ace  or  comma  between  the  data. 
AU  files  with  die  fine  extenraon  follow  die  same  format 
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ohmon.msc:  Block  1,  some  miscellaneous  data. 


Card  1  (Bank  rejection  options) 

Example: 

inst 


Variable 

name 

Type  and 
length 

Column 

number 

Description 

rejp 

a4 

Bank  rejection  parameter;  for  instantaneous  rejection,  = 
"inst;"  if = "full,"  then  the  oil  deposited  on  the  shoreline  dur¬ 
ing  each  time  step  will  be  rejected  completely  at  the  end  of 
the  time  step;  for  other  rejection  rates,  =  "part" 

Card  2  (Output  options) 
Example: 

0 

Variable 

name 

Type  and 
length 

Column 

number 

Description 

ipropt 

Integer 

— 

Two  options:  larger  than  1,  output  grid  lodty,  stream-tube 
velocity  to  data  files  for  plotting;  less  than  1,  no  output  to 
datafil^ 

Car4? 

Example: 

1.0 

Variable 

name 

Type  and 
Imgth 

Column 

number 

Description 

ajpha 

Real 

- 

Resurfacing  coefficient,  a;  if  all  particles  that  rise  to  the  sur¬ 
face  layer  stay  in  surface  layer,  =  1;  otherwise,  less  than  1 .0 

Card4 

Example: 

Le-5 

Variable 

name 

Typeand 

length 

Column 

number 

Description 

beta 

Real 

- 

Deposition  coefficient 

Example: 

0.001 

Vttriai^ 

name 

Typeand 

length 

Column 

number 

Description 

gama 

Real 

- 

Emulsification  coefficient  (1/s) 

Example: 

0.028 

VarkMe 

name 

Typeand 

kn^ 

Column 

number 

Description 

saxtiffi 

Real 

- 

Manning's  coefficient  of  riverbed 
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Card? 

Example; 

322 


Variable 

name 

Typeand 

length 

Column 

number 

Description 

gravac 

Real 

- 

Gravitational  acceleration  (ft/s^) 

Cards 

Example; 

0.01 

Variable 

name 

Typeand 

length 

Column 

number 

Description 

vbuoy 

Real 

- 

Buoyant  velocity  (ft/s) 

Card  9 

Example: 

1.1 

Variable 

name 

Typeand 

length 

Column 

number 

Description 

alpc 

Real 

— 

Coefficient  for  converting  depth-averaged  velocity  to  surface 
velocity, 

Card  10 

Example; 

0.03 

Variable 

name 

Typeand 

length 

Column 

number 

Description 

alpw 

Real 

- 

Wind  effect  coefficient  on  water  surface  velocity,  a^, 

Card  11 
Example: 
05  05 

Variable 

name 

Typeand 

length 

Column 

number 

Description 

ratsur 

ratsus 

Real 

Real 

Fraction  of  oil  spilled  on  water  surface 

Fraction  of  oil  spilled  in  the  suspended  layer 

^12 

Example: 

0.3 

Variable 

name 

Typeand 

length 

Column 

number 

Description 

alevp 

Real 

— 

Coefficient  for  modifying  evaporation  due  to  emulsification, 
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Card  13 
Example: 
1000  200. 


Variable  Type  and  Column 

name  length  number  Description 


nLxgrid 

dx 

Integer 

Real 

— 

Total  x-grid  of  moving  grids 

Grid  size  (ft) 

Card  14 
Example: 

5 

Variable 

Type  and 

Column 

name 

length 

number 

Description 

nr 

Integer 

- 

Total  number  of  river  segments 

ohmon.msc.  Blcxic  2,  coordinate  transformation  information. 

Cardl 

Example: 

129. 

Variable 

Type  and 

Column 

name 

length 

number 

Description 

sO 

Real 

- 

River  mileage  of  the  dovkmstream  boundary  of  the  study 
reach 

CwI2 

Example: 

0.0 

Variable 

name 

Type  and 
length 

Column 

number 

0.0  348200.0  144000. 

Description 

xmin(i) 

Real 

— 

Minimum  x<oordinate  of  i^  river  segment 

ymin(i) 

Real 

- 

Minimum  y-coordinate  of  i*  river  segment 

xmax(i) 

Real 

- 

Maximum  x-coordinate  of  river  segment 

yn)ax(i) 

Real 

- 

Maximum  y-coordinate  of  river  segment 

Card  3 
Example: 
0.0 

Variable 

name 

Type  and 
Imgth 

Column 

number 

0.0  0.0 

Description 

x0(i) 

Real 

- 

x-origin  (coordinate)  of  the  i*  local  coordinate  system 

y0(i) 

Real 

- 

y-origin  (coordinate)  of  the  i*^  local  coordinate  system 

beta0(i) 

Real 

- 

Rotation  angle  of  the  i*^  local  coordinate  x  positive  with  glo¬ 

bal  X  positive;  for  clockwise  and  "+"  for  counterclock¬ 

wise 
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Card  4 

Example; 

1741  46. 

Variable  Type  and  Column 

name  length  number  Description 

ngrid(i)  Integer  -  Maximijm  .t-grid  number  in  the  river  segment 

sbm(i)  Real  -  River  mileage  at  the  end  of  the ;  *  river  segment 

Note:  Repeat  cards  3  and  4  times  as  many  times  as  the  number  of  river  segments. 

This  is  the  end  of  data  file  ohmon.msc. 

Input  adjustments 

For  a  river  that  already  has  the  necessary  input  files,  very  little  modification  in  input  is  needed  to  run 
the  model  with  different  field  conditions,  such  as  new  river  discharge,  different  oil  properties,  new  spill 
location,  etc.  Input  data  that  are  most  likely  to  require  modification  are  cited  below  with  some  guide¬ 
lines  and  suggested  values.  Readers  should  refer  to  the  preceding  section  for  the  format  of  the  data. 

Velocity 

Interpolating  stream-tube  velocities  to  boxes  depends  on  the  variable  kintm  {phmon  .geo,  Block  1,  Card 
1)  and  ibend(i)  {ohmon.geo,  Block  2,  Card  1).  As  kintm  increases,  more  velocity  interpolations  are  made 
between  cross  sections.  If  the  cross  sections  are  spaced  far  apart,  kintm  should  be  in  the  range  of  5-8  or 
more;  otherwise  kintm  can  be  5  or  less.  The  prob¬ 
lem  with  using  too  small  a  value  for  kintm  is  that 
more  grid  boxes  will  end  up  without  assigned 
velocities  after  the  interpolation  stage.  They  will 
be  assigned  a  velocity  based  on  their  neighbor¬ 
ing  boxes.  This  is  less  accurate  than  assigning  a 
velocity  based  on  interpolatioit  One  should  try 
to  assign  as  many  boxes  as  possible  in  the  inter¬ 
polation  stage.  On  the  other  hand,  if  the  value  of 
kintm  is  too  large,  it  will  take  more  computation 
time.  Besides  the  value  of  velocities,  appropriate 
simulation  of  the  direction  of  velocity  is  also  im¬ 
portant.  Ibend(i)  supplies  two  options  for  deter¬ 
mining  velocity  directions.  When  ibendfi)  =  0,  the 
velocity  will  be  along  the  line  connecting  two  cor¬ 
responding  stream-tube  centers  in  neighboring 
cross  sections.  This  usually  gives  a  fairly  accurate 
flow  direction.  However,  for  cross  sections  in  a 
river  bend,  as  diown  in  Figure  26,  ibend(i)  should 
be  1.  This  naakes  the  velocity  direction  perpen¬ 
dicular  to  the  cross  section. 

Another  factor  that  affects  assigning  velocities 
to  grid  boxes  i  the  number  of  stream  tubes  se¬ 
lected.  The  choice  of  the  number  of  stream  tubes 
(nstube,  ohmon.geo,  Block  2,  Card  1)  will  affect  the 
ccanputer  time  and  the  accuracy.  For  a  river  with 
highly  irregular  cross  sectioris  or  when  a  high 
degree  of  accuracy  is  required,  more  stream  tubes 
i^Knild  be  used.  Figure  26.  Simulation  of  flow  direction  in  stream  tubes. 
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Stage-discharge  data 

Stage  and  discharges  (ohtnon.flo.  Block  1,  Card  2)  in  river  branches  are  unlikely  to  remain  constant 
over  time.  Therefore,  this  information  will  require  updates  as  the  need  arises.  The  one-dimensional 
unsteady  flow  model  can  be  used  to  obtain  this  information. 

Oil  spill  parameters 

The  oil  spill  parameters  (ohmon.spl,  Block  1,  Cards  2, 3, 4, 5)  will  require  the  most  modification  from 
spill  to  spill.  Each  spiU  will  have  its  own  simulation  duration,  simulation  time  step,  printing  options, 
number  of  particles,  volume,  physical  properties  and  location,  so  it  will  be  necessary  to  adjust  these 
data  each  time  the  simulation  is  done  for  a  new  spire  The  suggested  guidelines  for  these  parameters 
follow. 

The  initial  total  particle  number,  ntotal,  is  about  500;  a  maximum  of  1000  is  possible  in  the  current 
versioa  During  the  simulation  the  number  may  increase  to  several  thousand  because  of  the  exchange 
between  the  two  layers.  A  larger  particle  number  will  give  a  better  result  at  the  expense  of  a  longer 
executicHi  time.  The  time  step  suggested  for  spildt  is  900  s  (15  minutes). 

For  other  spill  parameters  the  following  values  may  be  used; 
spgoil  =  0.7-0.98 
aniu  =  1.411  x  lO*^  ft2/ s  (at  50°F) 
sigma  =  2.06  x  10-5  (]b/ft) 
ak2i  =  1.14 
ak2v  =  0.98 
ak2t  =  1.6 
akcl0  =  1.39 
akc20  =  1.39 
akc30  =  1.43 

vmnui  =  0.7063  x  10-2  (ft^/mol) 
soluni  =  0.1873  x  10-2  ab/ft3). 

Wind  data 

The  computer  model  reads  the  wind  data  iohmon.spl,  Block  2,  Card  1)  at  each  time  step.  This  allows 
for  changing  the  wind  speed  and  direction  in  the  simulation.  A  correction  is  needed  because  the  wind 
speed  on  land  and  that  on  the  river  are  different.  If  no  better  information  is  available,  a  value  of  0.2  may 
be  used  for  vwc  {ohmon.spl.  Block  2,  Card  1). 

Sample  iiq>ut  data  files 

All  file  names  to  be  opened  for  input  and  output  will  be  in  flie  file  ohmon.fnm.  Sample  input  data  files 
included  in  fliis  section  are: 

1.  (^tmofugeo 

River  geometry,  cross  sections  and  branch  geometry  of  the  C^o-Monongahela  rivers. 

2.  chnton.igr 

Boundary  grid  data  for  the  Ohio-Moncmgahela  rivers. 

3.  ohmotuice 

Ice  parameters,  ice  regions  and  ice  thickness  in  the  Ohio-Monongahela  rivers.  This  file  is  needed 
even  when  foe  model  is  run  for  open-water  conditions.  Samples  will  be  given  for  the  ice-covered 
condifion  and  foe  open^watar  conditiorL 

4.  ohmmflo 

Water  level  and  discharge  at  nodes  from  the  unst^dy  flow  model  for  the  Obio-Monongahela 
rivers. 

5.  ohtmmbnd 

Hall-life  index  for  sluneline  grids  of  foe  Cfoio-Moncmgahela  rivers. 
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6.  ohmonspl 

Oil  parameters,  spin  Icxration,  wind  condition  and  ambient  air  temperature. 

7.  ohmon.ndc 

Information  for  converting  unsteady  flow  computation  results  to  what  are  needed  in  oil  spill  model 
flow  computation  for  the  Ohio-Monongahela  rivers. 

8.  ohmon.msc 

River  division  data,  exchange  parameters  for  the  two  layers  and  information  for  coordinating  trans¬ 
formation  data  for  the  Ohio-Monongahela  rivers. 

ohmon. fnmThu  Jan  24  08:59:15  1991  1 


ohmon . f nm; 


lohmon . geo 
2oilprt.out 
7ohmon.flo 
Sohmon.bnd 
9ohinon .  ndc 
1 lohmon. sur 
13ohinon .  sus 
12ohmon . spl 
14ohinon.  ice 
20ohmon.msc 
22concen.dat 
29ohmon. Igr 
28vel .pos 
SOdep.pos 


ohmon . geo : 


24  16 

12  22  32  44  45  64 

89  104  108 

region  1 

121  131  138  139 

141 

150  165 

177 

194 

204  208 

210  218  228  233 

1  16840.66  7741.73 

-0.44720 

0.89443 

11 

11 

2 

0 

0.0 

0 

2  18925.39  7018.47 

-0.44720 

0.89443 

11 

11 

1 

0 

0.0 

0 

3  20144.21  8923.44 

-0.56714 

0.82362 

11 

11 

2 

0 

0.0 

0 

1  10 


0.0  660.0 

540.0 

660.0 

540.0 

596.0 

600.0 

596.0 

650.0 

596.0 

700.0596.0 

800.0 

596.0 

900.0 

596.0 

900.0 

660.0 

1274.0 

660.0 

2  20 

0.0  665.0 

1200.0 

665.0 

1280.0 

660.0 

1380.0 

625.0 

1659.0 

625.0 

1800.0616.0 

1900.0 

594.0 

2200.0 

594.0 

2420.0 

580.0 

2770.0 

581.0 

2940.0583.0 

3080.0 

590.0 

3140.0 

590.0 

3240.0 

606.0 

3300.0 

606.0 

3380.0633.0 

3830.0 

635.0 

3870.0 

640.0 

4030.0 

645.0 

4500.0 

665.0 

0 

23  5  region  2 
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ohmon.igr  Thu  Jan  24  09:01:241991  1 
ohmon . igr 


1 

15 

19 

2 

15 

19 

3 

15 

19 

4 

15 

19 

5 

14 

19 

6 

14 

19 

7 

14 

19 

8 

14 

19 

9 

14 

19 

10 

14 

19 

11 

14 

19 

12 

14 

19 

13 

14 

19 

14 

14 

19 

15 

14 

19 

16 

14 

20 

17 

14 

20 

18 

14 

20 

19 

14 

20 

20 

14 

20 

21 

14 

21 

22 

15 

22 

23 

15 

22 

24 

15 

23 

25 

15 

23 

26 

16 

23 

ohmon 

-ice: 

0.035 

0.84 

0 

0.035 

0.84 

0 

0.03S 

0.84 

0 

0.035 

0.84 

0 

0.035 

0.84 

0 

0.035 

0.84 

0 

0.035 

0.84 

0 

0.035 

0.84 

0 

0.035 

0.84 

0 

0.035 

0.84 

0 

0-035 

0.84 
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0 


0.035  0.84 


0 

0.035  0.84 

0 

ohmon .  floThu  Jan  24  09:03:01  1991  1 

ohmon. flo: 

6.0  flow  time  step 

branch  qu  qd  hu  hd  region  1  1/  1/88  hour  =  0 

1  52810.  52810.  625.19  625.00 

2  51761.  52810.  625.25  625.19 

3  50721.  51761.  625.37  625.25 

23  48876.  48878.  669.09  668.89 

24  48877.  48876.  669.30  669.09 

1 

2  open 

branch  qu  qd  hu  hd  region  2  1/  1/88  hour  =  0 

1  48877.  48876.  669.30  669.09 

2  48860.  48877.  669.45  669.30 

3  48850.  48860.  670.50  669.45 

4  48460.  48850,  670.  87  670.50 


ohmon. bnd: 

1  1  81  1 
2  1  81  2 

1  82  84  2 

2  82  84  2 

3  82  84  2 

4  82  84  2 

1  85  34  2 

2  85  34  2 

1  343  358  2 

2  343  358  2 


3  343 

358 

2 

4  343 

358 

2 

1  359 

530 

1 

2  359 

530 

1 

1  531 

532 

2 

2  531 

532 

2 

3  531 

532 

2 

4  531 

532 

2 

1  533  536  1 

2  533  536  1 

1  537  538  2 

2  537  538  2 

3  537  536  2 

4  537  538  2 

1  539  740  1 

2  539  740  1 
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ohmon . spl Thu  Jan  24  09:04:27  1991  1 

ohmon . spl : 


Ohio  Monongahela 
No.  2  Fuel  oil 

240.0  288  1  1  1  03600.  -1.0  totime, 

400  705000.0  300.  84  1.411E-5  7.55E-4  .14  1.1  1.6  1.4 


342650 

.-97800. . 

7063E-02 .1873E-02 

0.2 

11.80 

210.50 

24.00 

0.2 

4.21 

100.50 

20.00 

0.2 

0.00 

0.00 

17.50 

0.2 

2.53 

40.50 

15.50 

0.2 

2.53 

40.50 

14.50 

0.2 

0.00 

00.00 

18.00 

0.2 

0.00 

00.00 

25.00 

0.2 

3.37 

60.50 

29.50 

0.2 

3.37 

60.50 

30.50 

0.2 

2.53 

90.00 

29.50 

0.2 

7.59 

200.50 

28.50 

0,2 

10.96 

230.00 

28,00 

0.2 

14.33 

230. SO 

27.50 

0.2 

28.66 

240.50 

27.00 

0,2 

36.25 

250.50 

27.50 

0.2 

37.93 

280.00 

25.50 

0.2 

33.72 

280. SO 

19.50 

0.2 

26.13 

270,50 

13.50 

0.2 

25.29 

280.00 

9.00 

ohmon.] 

ndc: 

Region  l;part  1 

1 

12  1 

2  3 

4 

5 

2 

10  13 

14  15 

16 

17 

3 

10  23 

24  25 

26 

27 

4 

12  33 

34  35 

36 

37 

5 

1  45 

22 

8  211 

212  213 

214 

215 

23 

10  219 

220  221 

222 

223 

24 

6  229 

230  231 

232 

233 

Next 

we  have 

located  x-section 

1 

2  2.11 

1  0.09 

2 

2  2.00 

1  0.20 

231 

25 

2.50 

24  1 

.00 

232 

25 

2.04 

24  1 

.46 

233 

25 

1.50 

24  2 

.00 

234 

25 

1.02 

24  2 

.48 

Region  2:  part  1 

1 

8  1 

2  3 

4 

5 

7.88  465.0 


6 

7 

8 

9 

10 

11  12 

18 

19 

20 

21 

22 

28 

29 

30 

31 

32 

38 

39 

40 

41 

42 

43  44 

216 

217 

218 

224 

225 

226  227 

228 

234 


6  7  8 


ievery, 
1.4  1.4 
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ohmon.msc  Sun  Jan  20  16:42:57 


1991  1 


part 

0 

1. 

l.e-5 

0.00001 

0.028 

32.2 

0.01 

1.1 

0.03 

0.8 

0.2 

1000 

5 

129. 

0. 

0. 

1741 

348200. 

348200. 

2621 

378400. 

378400. 

2861 

315800. 

392400. 

3244 

52000. 

315800. 

4563 

2288 


full  for  rejrat  =  1  , inst  for  inst .  rej.,  part  for  calculated; 
=>ipropt,  >1  output  vel.  and  depth  in  grid,  opposite  in  inverse. 
=alpha,  resurfacing 

=beta,  deposition 

=gainma,  emulsification 


=ainngn 

=gravac 

=vbuoy 

=alpc 

=alpw 

0.2  =ratsur,  ratsus 
=alevp,  evap.  modified  coeff. 

200.  =ntxgrid,  dx  the  x-grid  number  of  working  length 

=nr 


0. 

348200. 

0. 

0. 

46. 

-36800. 

432000. 

139200. 

-90. 

0.0 

-84800. 

400000. 

-36800. 

-90. 

-11.5 

-132000. 

392400. 

-84800. 

-180. 

-40.5 

-220000. 

315800. 

-104600. 

-180. 

-132. 

3288 

sO 

144000.  xmin,  ymin,  .xina.x,yinax 

xO, y0,beta0 
ngrid(l),  sbrn(l) 
139200.  xmin,  ymin,xiT)ax,yinax 

x0,y0,beta0 
ngrid(2),  sbrn(2) 
-36800.  xinin,ymin,xmax,ymax 

x0,y0,beta0 
ngrid(3),  sbrn(3) 
-84800.  xmin,  ymin,xmax,yirax 

x0,y0,beta0 
ngrid ( 4 ) ,  sbrn ( 4 ) 
-104600.  xmin, ymin, xmax,ymax 

xO, y0,beta0 
ngrid (5),  sbrn (5) 


Additional  input 

Besides  the  input  data  files  mentioned  above,  there  are  two  additional  sets  of  data  to  be  input  on  the 
soeen  after  the  start  of  the  simulation. 

Fiver  combinations 

The  model  is  designed  to  simulate  spills  ttiat  occur  either  in  the  Monongahela  River  or  in  the  Allegh¬ 
eny  River.  After  the  user  starts  running  the  oil  spill  model,  the  screen  will  show; 


*  ★★★★★★*AA*******<t**;*r****A******AA**** 

*  select  river  combination:  * 

*  for  Ohio-Monongahela  River  input:  1  * 

*  for  Ohio-Allegheny  River  input;  2  * 


Choose  1  for  an  oil  spill  in  the  Monongahela  River  or  2  for  a  spill  in  the  Allegheny  River.  If  the  spill  is 
in  the  Ohio  River,  dhoose  eiftter  1  or  2.  However,  2  will  make  take  less  computing  time,  since  the  data 
files  on  ftte  geometry  and  shoreline  boundary  conditions  are  smaller. 


57 


Continuation  runs 

The  model  can  be  used  either  for  a  simulation  from  the  beginnirig  of  a  spill  or  for  continuing  a  previ¬ 
ously  stopped  simulation.  The  user  will  be  required  to  choose  from  these  two  options  after  choosing 
the  river  combination.  The  screen  wUi  show  the  following; 


*  The  model  can  be  run  from  t  =  0  or  * 

*  restarted  after  a  stop  * 

*  Enter  0  or  1  and  hit  return  * 

*  0  for  a  run  starting  from  t  =  0  * 

*  1  to  start  from  a  previously  stopped  point  * 

*  - - -  ★ 


MODEL  OUTPUTS 

Output  from  the  model  can  be  directed  to  different  devices,  such  as  a  console,  a  printer  or  files. 

Velocity  and  depth  distributions 

Flow  velocity  and  depth,  both  along  the  stream  tubes  and  in  the  grids,  are  calculated  in  subroutine 
xxldep.  If  needed  the  user  can  choose  the  printout  option  (ipropt>  1,  ohmon.int,  Block  1,  Card  2)  to  save 
all  simtilation  results  into  corresponding  files. 

•  strm.vel:  Save  coordinate  (x,y)  and  velocity  (v^,  Vy)  at  stream-tube  locations. 

•  strm.dep:  Save  coordinate  (x,y)  and  depth  (ft)  at  stream-tube  locations. 

•  grid.vel:  Save  coordinate  {x,y)  and  velocity  {v^,  Vy)  at  grid  box  locatioi\s. 

•  grid.dep:  Save  coordinate  {x,y)  and  depth  (ft)  at  grid  box  locations. 

All  of  these  data  are  used  in  advection  and  mixing  calculations.  They  cannot  be  saved  in  the  rniao- 
comjjuter  version  because  of  the  limitation  on  the  number  of  files  to  be  opened.  The  arrangement  of 
data  files  is  shown  as  follows: 

Data  card  of  velocity  for  grids  and  stream  tubes: 

X  y  Vy 

Data  card  of  depffi  for  grids  and  stream  tubes: 

X  y  h 

where  coordinates  and  depth  are  in  feet  and  velocity  is  in  feet  per  second. 

The  cards  for  velocity  and  depth  are  repeated  as  nrany  times  as  the  distribution  points.  In  this  mod¬ 
el,  only  the  velocities  and  depths  within  the  moving  grid  at  each  flow  time  step  will  be  saved  into  files 
to  save  memory  and  computer  time. 

Simulation  results 

Oil  spill  simulation  results  will  be  output  to  files  on  the  selected  output  ft’equency  (ieveiy ,  ohmon.spl. 
Block  1,  Gird  2).  These  results  include  the  following  files: 

•oUprt.out:  Spill  site  and  oil  parameters.  A  summary  of  distributions  of  the  oil,  including  updated 
total  volume  (for  checking  purposes),  volvime  in  the  surface  layer,  volume  in  the  suspended  layer, 
total  volume  evaporated  and  evaporated  in  the  current  time  step,  total  volume  dissolved  and 
volume  dissolved  during  this  time  step,  volume  on  land,  volume  deposited  onto  the  riverbed,  and 
vdume  emubified  and  resurfoced. 

•ohmonmr:  Locations  and  volumes  of  surface  oil  particles,  including  updated  total  particle  number 
and  foe  center  of  slick. 
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ohmon.sur:  Block  1 


Card  1 
Example: 

Ohio  Monongahela  Surface  No.  2  Fuel  Oil 


Variable 

name 

Type  and 
length 

Column 

number 

Description 

text 

a20 

11 

River  name 

sur 

al4 

31 

Indicate  results  of  surface  slick 

ftiellp 

a20 

45 

Name  of  spilled  oil 

Card  2 

Example: 

881.25 

Variable 

Type  and 

Column 

name 

length 

number 

Description 

vave 

- 

Average  volume  of  particles  both  on  surface  and  in  suspen¬ 
sion  at  the  begirming  for  graphics  purpose  (gal.) 

Cards 

Example: 

400  705000.  900. 

0.81  0.141E-04  0.755E-03  0.706E-02  0.187E-02  O.OOOE+00 

564000.000 

141000.000 

342650 

-97800 

Variable 

Type  and 

Column 

name 

Imgth 

number 

Description 

ntotal 

Integer 

Total  number  of  particles  defined  in  each  layer  (e.g.  400) 

spvol 

Real 

- 

Total  volume  of  oil  spilled 

spgoil 

Real 

- 

Specific  gravity  of  oil 

aitiu 

Real 

- 

Kinematic  viscosity  of  water  (ft^/s) 

sigma 

Real 

- 

Surface  torsion  of  oil  (Ib/ ft) 

vmuni 

- 

Molar  voliime  of  oil  (ft^/ mol) 

soluni 

- 

Solubility  of  fresh  oil  (Ib^/ ft^) 

amuni 

- 

Absolute  viscosity  of  oil  (lb/ ft^) 

Total  volume  of  oil  on  water  surface  at  the  beginning  of  spill 
(gal) 

Total  volume  of  oil  in  suspensicm  at  the  beginning  of  spill 
(gal.) 

ohmonsur:  Block  2 


Cardl 
Example: 
342650  -97800 


Variable 

name 

Type  and 
length 

Column 

number 

Description 

ipartx(l) 

i8 

_ 

x-coordinate  of  the  spill  site  (ft) 

iparty(l) 

i8 

9 

y<oordinate  of  the  spill  site  (ft) 

Card  2 

Example: 

2300  391145  -77512 

24.0000  - 

■1.6  -0.5 

Variable 

name 

Type  and 
length 

Column 

number 

Description 

nwtl 

Integer 

Total  number  of  surface  oil  particles  at  the  present  time 

ipartx(l) 

i8 

- 

x<oordinate  of  surface  slick  center  (ft) 

iparty(l) 

i8 

9 

y<oordinate  of  surface  slick  center  (ft) 

ttt 

Real 

- 

Completed  simulation  time  (hr) 

vwx 

Real 

- 

Present  x  component  of  wind  speed  (mph) 

vwy 

Real 

- 

Present  y  component  of  wind  speed  (mph) 

Example: 

391163  -85810  53835  391160  -85399  150.81  391185  -84711  163.11391370 


-85965  132.45 

This  card  has  four  sets  of  values  [x,y,  vopsur(i)]  to  indicate  the  coordinates  and  volume  of  each  particle; 
X  and  y  in  are  in  feet,  and  the  volume  of  the  surface  particle  vopsur(i)  is  in  gallons.  This  card  will 
repeat  for  all  particles. 

Block  2  will  be  repeated  as  many  times  as  desired  output  times. 

•okmonsus:  Locations  and  volumes  of  suspended  oil  particles,  including  the  updated  total  particle 
numbor  and  the  spill  center.  The  output  format  is  the  same  as  ohmon.sur. 

•concenMat:  One-dimensional  cross-section  averaged  concentration  distribution  (including  the  sur¬ 
face  layer  and  tihe  suspended  layer),  indicating  the  volume  concentration  at  different  river  miles. 

coruxnMat:  Block  1 

Cardl 

Exan^le: 

&  time,  iO,  disOs  24.0000  ^  33.9545 

where  &  is  for  a  gr^fdtks  purpose;  time  is  the  total  lengfo  of  the  completed  simulation  (hr);  iO  is  the 
inhial  x-giid  number  of  the  moving  grids;  disO  is  foe  river  mile  of  foe  downstream  edge  of  the  moving 
grids  (mBe). 
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Card  2 
Example: 

2.827  0.5278E-02 


Vamble 

Type  and 

Column 

name 

length 

number 

Description 

dis 

Real 

- 

River  mile  at  which  cross-section-averaged  volume  concen¬ 
tration  is  calculated 

conc2 

Real 

- 

Volume  concentration  at  river  nule  distance  (10"3) 

Block  1  will  be  repeated  as  many  times  as  the  designed  output  number. 

•  concensur.  One-dimensional  surface  concentration  distribution,  indicating  the  volume  concentration 
at  different  river  miles  along  the  river.  This  file  will  not  be  available  in  the  output  for  a  micro¬ 
computer.  The  format  of  this  file  is  the  same  as  concen.dat. 

The  graphics  output  that  corresponds  to  these  outputs  can  be  obtained  by  using  the  graphics  pro¬ 
gram. 

Sample  output  files 

oilprt.out  Wed  Jan  30  11:53:42  19911 

oilprt.out: 


1  Ohio  Monongahela 

*  continuous  spill  * 

*at  * 

*342650. ,-97800.  * 

*  for  60.  min.  * 

simulation  period  “  240.0  hrs 


characteristics  of  spill 

no.  of  pairticles  : 

oil  spilled  : 

dt  for  spill  simulation  : 
specific  gravity  of  oil  : 
kinematic  visco.  of  water: 
surface  tension  : 

spreading  coefficients 
k2i  k2v  k2t 

1.14  1.10  1.60 

molar  volume  ; 

solubility  of  fresh  oil  : 
viscosity  of  oil  : 

Manning's  roughness  of  ice 


400 

705000.  gals  of  No.  2  Fuel  Oil 
300.  secs. 

0.84  (api  index  =  37.0) 
0.1411E-04  sq  ft/sec 
0.7550E-03  Ibs/ft 

clO  c20  c30 

1.40  1.40  1.40 

0.7063E-02  cu  ft/mol 
0.1873E-02  Ibs/cu  ft 
0.00  gms/cm-sec 
0.000 


api  option  is  not  selected— evap .  constants  are  c=  7.88 t0=  465.0 

surface  diffusion— default  formulation  is  used 

time,  iO,  lend,  nhitb= 

8.00  2064  3064  1346 

slick  condition  at  the  end  of  this  time  step 

total  volume  (cu.  ft)  =  94244.5 

total  volume  on  land  (cu.  ft)  =  44775. 

total  on  surface  =  46411.  total  in  suspension  =  714.53 

evaporation  (cu.  ft)  this  step  =  18.694total  evap. (cu.  ft)  =  2222.5 

dissolved  this  step  =  0.58905  total  =  29.087 

total  resurfacing  (cu.  ft)  =  .12147E+06 

total  deposition  (cu.  ft)  =  121.39 

total  emulsification  (cu.  ft)  =  13582. 

total  out  of  range  on  surface  (cu.  ft  =  .66263E-02 

total  out  of  range  in  suspended  (cu.  ft)  =  .54000E-02 

ooncen.dat  Thu  Jan  24  09:12:37  1991  1 

concen . dat : 

6  time,  iO,  dis0=  24.0000  2059  31.5455 

0.4180  0.5278E-02 
0.3801  0.1285E-03 
0.3378  0.1198E-01 
0.2152  0.692eE-04 
0.1616  0.4378E-01 
0.1193  0.1687E-01 
0.7693E-01  0.3473E-01 
0,2336E-01  0.1634 

-0.1899E-01  0.1171 

-0.7256E-01  0.3773 

-0.1261  0.1599 

-0.1944  0.9520E-01 
-0.2627  0.7579E-01 
-0.8132  0.7076E-01 
-0.85C6  0.7076E-01 
-0.9576  1.052 

-1.060  0.9058E-01 
-1.128  0.1246 

-1.213  0.2245 

-1.281  3.160 

-1.869  1.859 

-1.937  2.506 

-2.006  1.563 

-2.059  2.579 

-2.113  4.183 

-2.166  0.9436E-01 
-2.274  0.4299E-03 
-2.327  2.119 

-2.381  0.1234E-04 
-2.419  1.403 

-2.472  0.1794E-01 
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-2.510  0.6395 
-2.564  0.166" 
-2.602  0.3221 
-2.639  0.2070 
-2.677  0.1014 
-2.715  0.5485E-01 
-2.753  0.6229E-01 
-2.791  0.7658E-01 
-2.833  0.1138 
-2.871  0.4878E-01 
-2.909  0.1256E-01 
-2.947  0.2762E-01 
-3.000  0.5685E-01 


ohraon.surThu  Jan  24  09:10:02  1991 

ohmon. sur : 


1 


Ohio  Monongahela  Surface  No.  2  Fuel  Oil 


881.250 


400  705000.  300.  0.84  0 

564000.000  141000.00 


342650 

-97800 

3228 

383560 

-63631 

24.0000 

385603 

-55130 

127.61 

385022 

385453 

-55484 

90.32 

385005 

384204 

-56198 

145.88 

384541 

384681 

-55980 

115.19 

384548 

384164 

-56405 

92.04 

384390 

384212 

-56492 

115.06 

383893 

383831 

-56706 

102.75 

384108 

383958 

-56767 

121.98 

384161 

383848 

-57152 

91.06 

383700 

384393 

-56345 

132.88 

383628 

385530 

-55327 

90.76 

383568 

141E-04  0.75'E-030.706E-02 


0.3 

-1.7 

-55737 

151.35 

384943 

-55816 

-55790 

130.41 

384217 

-56667 

-S5888 

99.56 

385486 

-55416 

-56119 

218.39 

384121 

-56382 

-56438 

167.40 

384514 

-S6237 

-56521 

114  22 

384098 

-56850 

-56482 

108.71 

384059 

-56695 

-56739 

100.85 

383714 

-57187 

-57115 

100.12 

383734 

-56833 

-57093 

152.93 

383296 

-57726 

-57324 

304.25 

384432 

-56297 

0.187E-020.000E+00 


122.80 

385406 

-55514 

137.42 

167.51 

384269 

-56384 

102.60 

225.97 

384341 

-56130 

230.82 

91.33 

384184 

-56360 

162.97 

216.95 

384947 

-55902 

149.75 

92.86 

384082 

-56556 

170.92 

209.84 

383949 

-56903 

500.66 

114.55 

384053 

-56653 

131.79 

118.72 

384327 

-56450 

119.88 

301.19 

383217 

-57443 

321.61 

96.78 

383300 

-57414 

192.12 

ohmon, sus: 


Ohio  Monongahela  Suspended  No.  2  Fuel  Oil 
881.250 


400  705000.  300.  0.84  0.141E-04  0.755E-030.706E-02  0.187E-020.000E+00 

564000.000  141000.00 

342650  -97800 


324  383388 

-63317 

24.0000 

0.3 

-1.7 

381100  -60059 

107.54 

380927 

-62207 

103.10 

381289 

-62940 

118.05 

381496 

-63364 

88.44 

381612  -63445 

93.64 

383286 

-64942 

102.53 

3859C0 

-54920 

1.47 

385700 

-55100 

1.57 

385500  -55300 

1.23 

385500 

-55500 

4.27 

385100 

-55700 

6.67 

385300 

-55700 

0.57 

384500  -55900 

1,30 

384700 

-55900 

9.48 

384900 

-55900 

1.15 

384300 

-56100 

6.66 

384500  -56100 

9.24 

383900 

-56300 

0.08 

384100 

-56300 

11.60 

384300 

-56300 

5.42 

384500  -56300 

1.51 

383700 

-56500 

0.05 

383900 

-56500 

3.79 

384100 

-56500 

3.25 

384300  -56500 

1.67 

383500 

-56700 

0.06 

383700 

-56700 

5.42 

383900 

-56700 

21.05 

384100  -56700 

3.49 

384300 

-56700 

0.52 

383500 

-56900 

9.22 

383700 

-56900 

13.31 

383900  -56900 

5.42 

384100 

-56900 

0.58 

383300 

-57100 

20.31 

383500 

-57100 

13.78 

383700  -57100 

10.03 

383900 

-57100 

0,55 

383300 

-57300 

21.70 

383500 

-57300 

7.66 

383700  -57300 

3.99 

383100 

-57500 

25.84 

383300 

-57500 

18.08 

383500 

-57500 

1.74 
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CONCLUSION 


This  report  describes  the  model  formulation  and  implementation  of  a  two-layer  model,  ROSS2,  for 
simulating  oil  spills  in  rivers.  The  model  considers  the  oil  in  the  river  as  consisting  of  both  the  surface 
slick  and  suspended  oU  droplets  mixed  over  the  depth  of  the  flow.  The  oil  transformation  processes 
considered  in  the  model  include  ad\  ection,  mechanical  spreading,  turbulent  diffusion  and  mixing,  evap¬ 
oration,  dissolution,  shoreline  deposition  and  sinking. 

The  model  can  be  used  for  instantaneous  or  continuous  spills,  either  on  or  under  the  water  surface 
in  rivers  with  or  without  an  ice  cover.  Although  it  is  developed  for  simulating  oil  spills,  the  model  can 
be  applied  to  spills  of  other  hazardous  materials. 

The  model  has  been  implemented  for  the  Ohio-Monongahela-Allegheny  river  system  and  the  up¬ 
per  St.  Lawrence  River.  A  case  study  for  the  spill  on  the  Monongahela  River  in  January  1988  is  present¬ 
ed  along  with  detailed  explanations  of  the  program  structure  and  its  input  and  output.  The  simulation 
results  compared  well  with  field  observations. 
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APPENDIX  A:  DEFINITION  OF  FORTRAN  VARIABLES 

Die  program  was  written  in  FORTRAN  77  and  was  tested  on  the  f77  compiler.  In  this  appendix, 
important  FORTRAN  variables  are  given.  For  more  detailed  instructions  on  the  usage  of  the  program, 
refer  to  the  main  text.  Descriptions  of  the  main  program,  subroutines  and  preparation  of  data  files  are 
also  given  in  the  text.  The  program  listing  is  provided  in  Appendix  B. 


Variable 

Symbol 

Type 

Description 

akclO 

reed 

Fay's  coefficient  for  gravity-inertial  phase  one-dimensional 
spreading 

akc20 

^20 

real 

Fay's  coefficient  for  gravity-viscous  phase  one-dimensional 
spreading 

akc30 

^30 

real 

Fay's  coefficient  for  surface  tension-  viscous  phase  one-dimen¬ 
sional  spreading 

ak2i 

^2i 

real 

Fay's  coefficient  for  gravity-inertial  phase  axisymmetric  spread- 

ak2t 

^2t 

real 

Ulg 

Fay's  coefficient  for  surface  tension-viscous  phase  axisymmetric 
spreading 

ak2v 

hv 

real 

Fay's  coefficient  for  gravity-viscous  phase  axisymmetric  spread- 

alevp 

se 

real 

Ulg 

Coefficient  for  modifying  evaporation  due  to  emulsification  (-  0.2- 
0.5) 

aJpc 

Sc 

real 

Weighting  factor  for  the  rivei  current  velocity  to  compute  the  sur¬ 
face  drift  velocity  (- 1.1) 

alpha 

a 

real 

Probability  factor  for  re-entrainment  of  oO  globules  from  the  sus¬ 
pended  layer  to  the  surface  layer 

alpw 

real 

Weighting  factor  for  the  wind  velocity  to  compute  the  surface  drift 
velocity  {-  0.03) 

amngn 

n 

real 

Manning's  roughness  coefficient  for  the  riverbed 

angle 

real 

Acute  angle  of  the  major  axis  of  the  slick  with  the  x-axis  (radians) 

anice 

«i 

real 

Manning's  roughness  coefficient  for  ice  at  the  ice/ water  interface 

amiuo 

real 

Absolute  viscosity  of  oil  (gm-cm/ s) 

amuni 

real 

Absolute  viscosity  of  oil  (Ibf-s/ft^) 

aniu 

IV 

real 

Kinematic  viscosity  of  water  (ft^/ s) 

api 

AH 

real 

American  Petroleum  Institute  index  for  oils 

apitem 

veal 

Temporary  storage  for  API  values 

b^ 

P 

real 

Fractional  empirical  factor  for  bottom  deposition  of  oil 

cevp 

C 

real 

Evaporation  coefficient  C 

cordlb(m) 

complex 

Complex  variable  giving  x  and  y  coordinate  locating  cross  sec¬ 
tion  on  reference  shore 

cordv(m,n) 

complex 

Coordinates  at  which  vstrm(m4t)  is  acting 

ddd 

Dy 

real 

Horizontal  diffusion  coefficient 

depth(n) 

h 

real 

Eteptii  of  the  river  at  the  n*  grid  box  center  (ft) 

dstatnfm^t) 

real 

Flow  dq7ffi  at  cross  section  and  «•*'  stream  tube  (ft) 

difx 

Dx 

real 

Diffiision  coefficient  in  the  x  direction 

dify 

Dy 

real 

Diffusion  coefficient  in  the  y  direction 

dvdt 

Ov/dt 

real 

Rate  of  drange  of  oil  volume 

dx 

real 

Widdi  of  the  grid  in  x  or  y  directicm  (since  square  grids  are  used, 
tihe  widths  are  equal)  (ft) 
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fevpl 

real 

Fraction  of  oil  evaporated  at  the  end  of  previous  time  step 

fevp2 

real 

Fraction  of  oil  evaporated  at  the  end  of  present  time  step 

fueltp 

character 

Type  of  oil 

fuU 

character 

Character  variable 

gama 

Y 

real 

Empirical  fraction  for  oil  exchange  from  surface  to  suspended 
layer 

gravac 

8 

real 

Gravitational  accelerahon  (ft/s^) 

hlife{17) 

X 

real 

Half-life  of  the  shoreline  for  oil  retenhon  (hr) 

ibend(n) 

integer 

Index  for  the  cross  section;  if  0,  flow  direction  will  be  along  the 
connection  line  of  corresponding  stream  tubes  in  two  neighbor¬ 
ing  cross  sections;  if  1,  flow  direction  (shore  or  island) 

idam(n) 

integer 

Particle  numbers  that  have  hit  locks  and  dams;  stored  in  the  ar¬ 
ray  where  n  is  from  1  to  ndam 

idum(n) 

integer 

Temporary  storage 

ievery 

every  otfier  time  step 

integer 

Frequency  of  obtaining  output;  e.g.  value  of  2  gives  output 

igrilMi) 

integer 

Water-side  }  grid  number  (y  direction)  of  the  lower  boundary  of 
the  river  for  the  grid  in  the  x  direction 

igriub(i) 

integer 

Water-side  /  grid  number  (y  direction)  of  the  upper  boundary  of 
the  river  for  the  i*  grid  in  ^e  x  direction 

igrlb(k4) 

integer 

Land-side  j  grid  number  (y  direction)  of  lower  boundary  of 
island  (if  any)  for  die  f*  grid  in  the  x  direction 

igrub(k4) 

integer 

Land-side  ;  grid  number  (y  direction)  of  the  upper  boundary  of 
k*  island  (if  any)  for  the  grid  in  the  x  direction 

ihitb(n) 

integer 

Particle  numbers  diat  have  hit  land;  stored  in  the  array  where  n 
isfroml  tonhitb 

imovin(n) 

integer 

Particle  numbers  that  are  in  the  water;  stored  in  the  array  where 
n  is  from  1  to  nmovin 

indpm 

integer 

Index  for  printout;  if  0,  no  printout;  if  1,  printout 

inicbdd 

integer 

Index  for  spreading;  if  0,  axisymmetric  spreading;  if  1,  one-di- 
mensicmal  spreading 

ipos 

integer 

Grid  box  number 

ipoll(n) 

integer 

Beginning  grid  box  number  of  the  ice  region 

ipos2(n) 

integer 

Ending  grid  box  number  of  the  n*  ice  region 

izibx(n) 

integer 

X  grid  number  of  the  n**'  grid  box  to  have  zero  velocity 

izri^n) 

integer 

y  grid  number  of  the  n*  grid  box  to  have  zero  velocity 

jduin(ll) 

integer 

Temporary  storage 

kintm 

integer 

Number  of  interpolations  between  two  velocity  points  of  two 
consecutive  cross  sections  in  tire  same  stream  tube 

kstsq(m) 

integer 

Last  cross-section  number  of  branch  m 

Ihnax 

integer 

Maximum  i*  grid  of  the  present  moving  grid 

Umin 

integer 

Minimum  i’l'  grid  of  the  present  moving  grid 

Imax 

integer 

Maximum  grid  of  the  present  oil  slick 

lam 

inte^ 

Minimum  grid  of  the  present  oil  slick 

nbmdt 

integer 

Number  of  brartches  used  in  oil  spill  model 

ndam 

integer 

Total  number  of  particles  diat  have  hit  locks  and  dams 

nfix«o(n) 

integer 

Next  cross  section  connecting  to  cross  section  in  question.  For  a 
divided  channel  around  island,  his  represents  he  first  cross  sec¬ 
tion  ccamected  to  in  he  lower  division  horn  he  main  channel 
cross  section 

ngiiidx 

int^go’ 

Total  riuiiiber  of  grids  in  tlffi  X  direction 

integer 

Total  number  of  particles  hat  have  hit  larvl 

integer 

Total  number  of  ice  regicns 

nicexl(n) 

integer 

X  grid  box  of  the  starting  point  of  ice  region  n 

mcex2(n) 

integer 

X  grid  box  of  the  ending  point  of  ice  region  n 

niceyl(n) 

integer 

y  grid  box  of  the  starting  point  of  ice  region  n 

nicey2(n) 

integer 

y  grid  box  of  the  ending  point  of  ice  region  n 

nmovin 

integer 

Total  number  of  moving  particles 

nseco(n) 

integer 

For  a  divided  channel  around  an  island,  this  represents  the  first 
cross  section  connected  to  in  the  upper  division  from  the  main 
channel  cross  section  (if  no  island,  =  0;  if  lower  division  complete 
and  returning  back  to  upper  division,  =  888;  if  both  divisions  are 
complete  and  resuming  main  channel,  =  999) 

nsisct(n) 

integer 

Number  of  sounding  depths  used  to  describe  the  channel  geom¬ 
etry  at  the  cross  section 

nstube 

integer 

Total  number  of  stream  tubes  at  the  cross  section 

ntpar 

integer 

Total  number  of  particles 

numcon(n) 

integer 

Condition  number  of  cross  section  n.  If  all  stream  tubes  continue 
to  next  cross  section  imdivided,  =  11;  if  stream  tubes  divide  into 
two  chaimels  from  main  charmel,  =  12;  if  stream  tubes  from  di¬ 
vided  channels  cormect  back  to  main  channel,  =  21 

nzrvb 

integer 

Total  number  of  boxes  to  assign  zero  velocities 

parsur(n) 

complex 

X  and  y  coordinate  location  of  tiie  particle  in  the  surface  layer 
(ft,  ft) 

parsus(n) 

complex 

X  aird  y  coordinate  location  of  the  particle  in  the  suspended 
layer  (ft,  ft) 

preldt 

real 

Time  step  for  releasing  new  particles  in  case  of  continuous  spill 
(in  one  spildt  several  preldt  occur)  (s) 

qsct(n) 

real 

Discharge  at  the  cross  section  (ft^/s) 

taciivts(n) 

real 

Distance  of  the  particle  from  the  spill  centroid  in  case  of  axi- 
symmetric  spreadingor  the  distance  of  the  n*  particle  from  the 
strip  centroid  in  case  of  one-dimensional  spreading 

sctang(n) 

teal 

Angle  the  cross  section  makes  with  positive  x  axis  (radiara) 

sigma  a 

real 

Surfrice  tension  (Ibf  /  ft) 

slickr(n) 

real 

Slick  radius  of  the  n*  pie  segment  (ft) 

solblt 

real 

Solubility  of  oil  (gm/m^) 

soluni 

real 

Solubility  of  oil  (Ib/ ft^) 

spaice 

real 

Surface  area  of  spill  under  ice  (ft^) 

sparea 

real 

Open  surface  area  of  spiQ  (ft^) 

spcea 

complex 

X  and  y  coordinates  of  the  spill  centroid  (ft,  ft) 

spcenO 

complex 

X  and  y  coordinates  of  ttie  spill  centroid  at  the  begiruung  of  simu- 
laticm  (ft,  ft) 

spgoil 

real 

Specific  gravity  of  oil 

real 

c5il  spill  time  step  (s) 

splrat 

real 

Rate  of  oil  spill  in  case  of  continuous  spill  (ft^/s) 

spitim 

real 

Duration  of  oil  spill  in  case  of  continuous  spill  (s) 

spvol 

real 

Total  volume  of  spilled  (U.S.  gaL) 

text(n) 

character 

Storage  array  for  character  string 

tk£(m;n) 

real 

Thickness  of  ice  at  the  m*  cross  section  an  sounding  depth 
locatkm  (ft) 

timet 

real 

Total  time  from  the  beginning  of  simulation  to  the  end  of  present 
time  step  (s) 

totune 

real 

Total  time  for  whidh  tiie  simulation  has  to  be  performed  (hr) 

tolevp 

real 

Total  volume  of  oil  evaporated  imtil  the  end  of  present  time  step 
(ft3) 

Oil  rejection  rate  horn  the  ^ore.  The  dimension  4  is  for  the  four 

typbfict(44) 

real 
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shores  (lower:  1,  upper:  2,  lower  island:  3,  upper  island:  4).  The 
dimension  i  represents  the  grid  in  the  .r  direction 

tOevp 

real 

Initial  boUing  point  of  oil  (K) 

tOuni 

real 

Initial  boiling  point  of  oil  (°R) 

vbuoy 

V],  real 

Buoyant  velocity  of  the  oil  globules  in  the  suspended  layer  (ft/s) 

vcar(n) 

complex 

X  and  y  components  of  velocity  in  the  grid  box  (ft/ s,  ft /s) 

vmimi 

real 

Molar  volume  of  oil  (m^/ mol) 

vmol 

real 

Molar  volume  of  oil  (m^/  mol) 

volevp 

real 

Volume  of  oil  evaporated  during  the  present  time  step  (ft^) 

volpie(n) 

real 

Volume  of  oil  in  the  n*  pie  multiplied  by  8  (ft^) 

volspl 

real 

Total  volume  of  oU  spilled  (fP) 

vopsvir(n) 

real 

Oil  particle  volume  of  surface  particle 

vopsus(n) 

real 

Oil  particle  volume  of  suspended  particle 

volspl 

real 

Imaginary  original  volume  of  oil  spilled  (m^) 

vpsuiO 

real 

Initial  particle  volume  in  surface  layer  (ft^) 

vtsurO 

real 

Initial  total  volume  in  surface  layer  (ft^) 

vtsusO 

real 

Initial  total  volume  in  suspended  layer  (ft^) 

vstnn(m^) 

complex 

X  and  y  components  of  the  stream  velocity  at  the  m**’  cross  section 
and  stream  ube  (ft/s,  ft/s) 

vwind 

complex 

X  and  y  components  of  wind  velocity  (ft/ s,  ft/s) 

vwmag 

real 

Magnitude  of  the  wind  velocity  (ft/s) 

vzero 

real 

Volirme  of  the  spill  (ft^) 

wsct(n) 

real 

Water  level  at  the  n*  cross  section  (ft) 

word 

character 

Character  variable 

yshift(n) 

real 

Distance  the  particle  is  shifted  m  the  negative  y  direction  if 
there  is  an  island  (ft) 

ywid(m^) 

real 

Distance  along  the  m*  aoss  section  from  the  reference  bank  to 
ttte  w*  soimding  depth  location  (ft) 

z(m^) 

real 

n*  sounding  depth  for  the  m*  cross  section  (ft) 

zd(n) 

real 

Reference  datum  for  section  from  which  the  sounding  depth 
is  evaluated(ft) 
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APPENDIX  B:  LIST  OF  PROGRAMS 


The  computer  program  is  written  in  FORTRAN  77  and  has  been  tested  with  several  Fortran 
compilers.  The  program  usage  is  discussed  in  the  main  text  and  Appendix  A.  Model  output  (results 
files)  are  described  in  the  text.  The  main  program,  ROSS2,  with  all  of  the  following  subroutines  will 
be  saved  on  a  floppy  disk  with  input  and  output  sample  data  for  the  users'  reference. 


ross2 

advsur 

advsus 

areacl 

check 

checkl 

chif 

disolu 

emdam 

emulfn 

lumper 

ohndc 

orient 

presur 

presus 

random 

reloca 

veldep 

shift 

sprdld 

bndloc 

bndrej 

calcij 

ckout 

con 

conchk 

evapor 

exchan 

input 

pievol 

ploter 

positn 

rearrg 

redepn 

rejloc 

sprdax 

strvol 

ohvad 
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